Logical modeling of lymphoid and myeloid cell specification and transdifferentiation

Abstract

Blood cells are derived from a common set of hematopoietic stem cells, which differentiate into more specific progenitors of the myeloid and lymphoid lineages, ultimately leading to differentiated cells. This developmental process is controlled by a complex regulatory network involving cytokines and their receptors, transcription factors, and chromatin remodelers. Using public data and data from our own molecular genetic experiments (quantitative PCR, Western blot, EMSA) or genome-wide assays (RNA-sequencing, ChIP-sequencing), we have assembled a comprehensive regulatory network encompassing the main transcription factors and signaling components involved in myeloid and lymphoid development. Focusing on B-cell and macrophage development, we defined a qualitative dynamical model recapitulating cytokine-induced differentiation of common progenitors, the effect of various reported gene knockdowns, and the reprogramming of pre-B cells into macrophages induced by the ectopic expression of specific transcription factors. The resulting network model can be used as a template for the integration of new hematopoietic differentiation and transdifferentiation data to foster our understanding of lymphoid/myeloid cell-fate decisions.

Keywords: cell fate; cell reprogramming; dynamical modeling; gene network; hematopoiesis.

Fig. 1.

(A) Schematic representation of hematopoietic cell specification. Genes in red are required for progression at the corresponding steps. C/EBPα-induced transdifferentiation is indicated by red arrows from B-lineage cells to macrophages. (B) Iterative modeling workflow. A model is first built based on the literature and is used to predict dynamical behaviors (cell phenotype, differentiation, reprogramming, and so forth). Predictions then are compared with experimental data; when the predictions and experimental data agree, further predictive simulations are performed; when they do not agree, further regulations are inferred from ChIP-seq data and are integrated into the model until simulations fully agree with data.

Fig. S1.

(A) Heatmap showing the regulations inferred from the literature and from ChIP-seq meta-analysis. Green cells denote activations, red cells denote repressions, blue cells denote context-dependent regulations, and gray cells denote regulations supported only by physical evidence. The color intensity denotes the level of confidence of the regulation. (B) Gene-expression analysis by qPCR for Foxo1 (Upper) and Spi1 (Lower) in macrophages (RAW cell line) before and after ectopic expression of Foxo1. Expression was normalized with respect to GAPDH and to wild-type RAW cells. (C) Expression profiles of B-cell genes during transdifferentiation (measured by Affymetrix microarrays) (29). (D) Signal from ChIP-seq data targeting different TFs and the histone marker H3K27ac in B cells, in B cells after C/EBPα induction (B+C/EBPα), in GMPs and in macrophages at Pax5, MEf2c, and IL7r loci. Black frames indicate B-cell enhancers presumably inactivated upon Cebpa binding after its induction. The vertical axes represent RPM (maximum: 2 RPM for Ebf1, 5 RPM for other TFs, 3 RPM for H3K27ac).

Fig. 2.

(A) Heatmap showing the regulations inferred from the literature and from ChIP-seq meta-analysis. (B) ChIP-seq signals and peaks (under signal) at the Spi1 locus. Black frames indicate known enhancers (24). The vertical axes represent reads per million (RPM) (maximum: 2 RPM for Ebf1 and Ikaros, 1.5 RPM for Foxo1, 1 RPM for Runx1 and Gfi1, 5 RPM for other TF). (C) ChIP-seq signals and peaks (under signal) at the Foxo1 locus. Black frames indicate B-cell enhancers in which C/EBPα binding is detected. The vertical axes represent RPM (maximum: 2 RPM for Ebf1, 5 RPM for other TFs, 3 RPM for H3K27ac). Note that Pax5 and Ikaros peaks are located downstream of the first exon and all other peaks are upstream of the TSS.

Fig. 3.

A regulatory graph depicting the interactions inferred from the literature and ChIP-seq meta-analyses. Nodes represent genes (except for CSF1r_act and Il7r_act, which represent the activated forms of cytokine receptors), and arrows denote regulatory interactions. Orange nodes represent factors expressed in macrophages, purple nodes represent factors expressed in progenitors, and blue nodes represent factors expressed in B-lineage cells. Ellipses represent Boolean components; the rectangle emphasizes the ternary component Spi1. Green and red edges correspond to activations and inhibitions, respectively. Gray edges denote the regulations predicted by the ChIP-seq meta-analysis, which were included in the model to increase consistency with expression data.

Fig. 4.

(A) Gene-expression values (microarrays) in lymphoid/myeloid progenitors (LLPP), B cells, and macrophages (Mac) (29). These values are relative to the highest expression value. (B) Context-dependent stable states computed for the model. A yellow cell denotes the inactivation of the corresponding component, a red cell represents maximal activation (1 for Boolean components, 2 for Spi1), and an orange cell represents an intermediate level (1) for Spi1.

Fig. S2.

(A) Context-dependent stable states obtained for our first literature-based model computed using GINsim software. The columns list the stable states. A yellow cell denotes the inactivation of the corresponding component (row), a red cell represents a maximal activation (1 for Boolean components, 2 for Spi1), and an orange cell represents an intermediate level (1) for Spi1. Note that the values of the input nodes (Csf1, Il7) are omitted. (B) Signal and peaks (under signal) from ChIP-seq data targeting different TFs in B cells and macrophages at the Cebpa locus. Black frames indicate loci of Foxo1 binding (Upper) or of myeloid factors binding (Lower). The vertical axes represent RPM (maximum 2 RPM for Ebf1 and Ikaros, 1.5 RPM for Foxo1, 5 RPM for other TFs). (C) Gene-expression analysis by qPCR for Foxo1 (Left) and Cebpa (Right) in macrophages (RAW cell line) before and after ectopic expression of Foxo1. Expression was normalized with respect to GAPDH and to wild-type RAW cells. (D) Context-dependent stable states of our model, as in A, after the consideration of Cebpa repression by Foxo1. (E) Signal and peaks (under signal) of ChIP-seq data targeting different B-cell factors in B cells. The vertical axes represent RPM (maximum: 2 RPM for Ebf1 and Ikaros, 1.5 RPM for Foxo1, 5 RPM for other TFs). (F) RNA-seq data (47) for hematopoietic cells corresponding to the stable states of our model. Values are normalized with respect to the highest expression value. (G) Gene expression for Egr1, Egr2, and Gfi1 as measured with microarrays (Left) (29) and RNA-seq (Right) (47). The y axis represents normalized average probe intensity for microarray and reads per kilobase of transcript per million reads mapped (RPKM) for RNA-seq. (H) Western blot for EGR2, EGR1, GFI1, and GAPDH proteins in pre-B and macrophage cell lines (C10 and RAW, respectively). Note that EGR1 and GFI1 were detected on the same blot, but pictures were taken with different exposure times, so their respective GAPDH control lanes are similar. (I) Sensitivity analysis of the proportion of lymphoid and myeloid cells after the differentiation of MPs in the absence of cytokine stimulation. Each bar represents a single simulation in which the up-regulation (Upper) or down-regulation (Lower) rate of a single variable was changed by more than six orders of magnitudes. The blue section of the bar represents the proportion of CLPs, and the orange sections represent the proportion of GMPs. Simulations with default rates equal to 1 (as in Fig. 5B) are highlighted with a gray frame. Gray horizontal lines show the proportions of CLPs and GMPs with default parameters.

Fig. 5.

(A) State transition graph generated by simulating the model starting from the unstable MP state in the absence of cytokine (Upper) and after the addition of CSF1 and IL7 (Lower Left and Lower Right). Nodes denote states, and arrows represent transitions between states. (B) Stochastic simulations showing the evolution over time, before and after cytokine exposition, of the fractions of cells expressing specific macrophage factors (Top), B-cell factors (Middle), and cell-type signatures (Bottom). The x and y axes represent time (in arbitrary units) and fractions of positive cells, respectively. (C) Hierarchical transition graph corresponding to the state transition graph in A. Nodes represent clusters of states, and arrows denote the possible transitions between the clusters. The labels associated with the edges highlight the crucial transitions involved in the decision between B-cell and macrophage specifications. (D–F) Schematic representations and stochastic simulations of the effects of Cebpa knockout (D), Pax5 knockout (E), or Spi1 knockout (F) on the differentiation of MPs, compared with the wild-type situation in A and B. In the cartoons, the wild-type stable states (cell types) and transitions that are lost in each mutant are displayed using light gray arrows and shading. MP, B cells, and macrophages are represented in purple, blue, and red, respectively.

Fig. S3.

(A) State transition graph obtained by simulating a permanent C/EBPα induction in B cells. Nodes represent states (i.e., vectors of values for all model components), and edges denote enabled transitions between states. The blue node corresponds to B cells, and the orange node corresponds to macrophages. (B) Expression profile for Cebpa, Cebpb, and Spi1 during transdifferentiation as measured by Affymetrix microarrays (29). (C) Stochastic simulations of B-cell transdifferentiation upon permanent Cebpa induction showing the temporal evolution of the fractions of cells expressing CD19 or Mac1 markers (Top), B-cell TFs (Middle), and macrophage TFs (Bottom). The x and y axes represent time (in arbitrary units) and fractions of positive cells, respectively. (D) Hierarchical transition graph obtained by simulating a C/EBPα transient induction in Spi1-KO B cells. The gray node on the right depicts the basin of attraction of the state in which all the components are inactivated (black node). Ebf1 is inactivated in all states of this basin of attraction and cannot be reactivated. (E) Hierarchical transition graph obtained by simulating a C/EBPα transient induction in B cells constitutively expressing Pax5.

Fig. 6.

(A, Upper) Hierarchical transition graph of the simulation of B-cell transdifferentiation upon transient C/EBPα expression, taking into account all possible C/EBPα pulse durations. Nodes represent clusters of states, and arcs correspond to transitions between these clusters. (Lower) Cebpa⁻ states (rows) of the basin of attraction of the macrophage stable state. (B) ChIP-seq signals and peaks (under signal) in B cells (in blue, time point 0 h) and after induction of C/EBPα (at 3, 12, and 24 h). The vertical axes represent RPM (maximum, 5 RPM). (C) Stochastic simulations of the fraction of cells expressing different B-cell factors (Top Panel) and cell-population signatures (Lower Three Panels) during transdifferentiation upon permanent (Upper Two Panels) or transient (Lower Two Panels) C/EBPα ectopic expression. The corresponding induction durations (in arbitrary units) are indicated by the black lines above each panel.

Fig. 7.

Table summarizing the impact of selected perturbations (knockin or knockout) on B-cell transdifferentiation into macrophages (Mac) upon either a permanent or a transient induction of C/EBPα. Orange boxes represent macrophages, blue boxes B cells, gray boxes all 0 stable states or cell death. Two-color boxes denote alternative outcomes (stable states).