Computational Modeling of a Transcriptional Switch Underlying B-Lymphocyte Lineage Commitment of Hematopoietic Multipotent Cells

Abstract

Despite progresses in identifying the cellular mechanisms at the basis of the differentiation of hematopoietic stem/progenitor cells, little is known about the regulatory circuitry at the basis of lineage commitment of hematopoietic multipotent progenitors. To address this issue, we propose a computational approach to give further insights in the comprehension of this genetic mechanism. Differently from T lymphopoiesis, however, there is at present no mathematical model describing lineage restriction of multipotent progenitors to early B-cell precursors. Here, we provide a first model—constructed on the basis of current experimental evidence from literature and of publicly available microarray datasets—of the genetic regulatory network driving the cellular fate determination at the stage of lymphoid lineage commitment, with particular regard to the multipotent-B-cell progenitor transition. By applying multistability analysis methods, we are able to assess the capability of the model to capture the experimentally observed switch-like commitment behavior. These methods allow us to confirm the central role of zinc finger protein 521 (ZNF521) in this process, that we had previously reported, and to identify a novel putative functional interaction for ZNF521, which is essential to realize such characteristic behavior. Moreover, using the devised model, we are able to rigorously analyze the mechanisms underpinning irreversibility of the physiological commitment step and to devise a possible reprogramming strategy, based on the combined modification of the expression of ZNF521 and EBF1.

About the Authors:

Luca Salerno

Contributed equally to this work with: Luca Salerno, Carlo Cosentino

Affiliation: Laboratory of Biomechatronics, Department of Experimental and Clinical Medicine, University Magna Græcia of Catanzaro, Catanzaro, Italy

Carlo Cosentino

Contributed equally to this work with: Luca Salerno, Carlo Cosentino

* E-mail: [email protected]

Affiliation: Laboratory of Biomechatronics, Department of Experimental and Clinical Medicine, University Magna Græcia of Catanzaro, Catanzaro, Italy

Giovanni Morrone

Affiliation: Laboratory of Molecular Haematopoiesis and Stem Cell Biology, Department of Experimental and Clinical Medicine, University Magna Græcia of Catanzaro, Catanzaro, Italy

Francesco Amato

Affiliation: Laboratory of Biomechatronics, Department of Experimental and Clinical Medicine, University Magna Græcia of Catanzaro, Catanzaro, Italy

Introduction

The differentiation of hematopoietic cell is governed by several control mechanisms, at different levels and in a highly hierarchical manner, from hematopoietic stem cells (HSCs) to lymphoid multipotent progenitors (LMPPs) that can further develop toward lineage restricted progenitors (LRPs). In these mechanisms, some transcription factors play a key role in the determination of the hematopoietic lineage. Observing these mechanisms from a systemic viewpoint, the following questions arise: (i) Which are the upstream regulators of the transcription factors that determine the lineage of differentiation of a certain progenitor or immature cell? (ii) Which are the transcription factors that constitute a minimal set for the activation of the differentiation process in hematopoietic cells? (iii) Is the differentiation an irreversible event (i.e., once the cell is committed to a certain path, it cannot change), or is it possible to observe phenomena of plasticity? While these and other issues, related to the regulation of hematopoietic development and differentiation, clearly require extensive and time-consuming experimental investigations, the exploitation of mathematical models and associated theoretical analysis tools can significantly boost the comprehension of the lineage specification mechanisms.

The control of HSCs differentiation towards mature hematopoietic cells is a sequential process, consisting of several lineage-determining steps. Each step requires a decision mechanism, which is assumed to be highly regulated by a combination of transcription factors, epigenetic events, and extrinsic regulator cytokines [1]. Prior to commitment, it has been observed that many genes involved in HSCs differentiation are expressed at intermediate or basal levels [2, 3]. This “transcriptional priming” at the stage of early progenitors may induce the rapid deployment of transcription factors to implement the subsequent commitment mechanisms. In hematopoiesis there exist several lineage branch points with identified key transcription factors and external signals [4–6]. For instance, a particularly well studied subnetwork is the one involving the genes PU.1 and GATA-1, which underlies the erythroid-myeloid lineage-determination step and has been proved to exhibit both a commitment switch and priming features [7–12]. In many cases, switch-like gene circuits have evolved to realize a sort of cellular built-in memory, which can lead to phenotypic diversity [13].

Several models have been proposed to describe a late stage of B-lymphocytes differentiation into terminal B-cells: sustained exposure to CD40L has been suggested to direct germinal center B-cells toward the memory B-cell compartment [14]; a bistable switch arising from the coupled double-negative feedback loops involving BCL-6, BLIMP-1 and PAX5 forms the basis of the B-cell to plasma cell differentiation program and its disruption by dioxin [15]; in [16], the authors developed a kinetic model to quantitatively characterize B-cell exit from the germinal center phase and terminal differentiation into plasma and memory B-cells.

Differently from those related to the terminal phases, the regulatory mechanisms underlying other lineage-determination steps are still poorly understood: for instance, to our knowledge, no model of the early B-lymphopoiesis commitment switch has been devised so far. Thus, the present work focuses on the analysis of this important mechanism. In particular, we first develop a dynamical model in order to clearly identify a regulatory network of genes/transcription factors involved in cellular fate determination of B-lymphoid lineage specification, based on the experimental evidence reported in literature. A special emphasis is given to the role of ZNF521, a key regulator and crucial antagonist of some genes involved in the B-lymphoid commitment regulatory network, as recently proved in experimental studies [17, 18]. Morrone and coworkers have suggested that ZNF521, in concert with other transcription factors, influences self-renewal and differentiation of primitive progenitors in hematopoiesis.

Exploiting the proposed mathematical model, we have studied the gene network regulating B-lymphoid commitment switch; in particular, we have investigated the role of several novel putative interactions by comparing the behavior of the alternative models against the available experimental data. Bifurcation analysis and numerical simulations have shown that the currently known interactions are not sufficient to generate the bistable behavior that distinguishes biological switches; moreover, the analysis of the extended models quantitatively proves that a novel interaction, that is the inhibition of ZNF521 expression by the B-cell lineage specific activator protein PAX5 [19], is required for the model to exhibit the commitment switch. Finally, the devised mathematical model is exploited to investigate possible strategies for reprogramming B-lymphoid LRPs into LMPPs.

Results

B-cell lineage commitment model and the role of ZNF521

Robust cell fate commitment is based on the activation and inhibition of the expression of certain genes. The two cytokine receptors, FLT3 and IL-7R, and the seven transcription factors PU.1, IKAROS, GFI1, E2A, EBF1, PAX5 and ZNF521 are critical for the development of B-cell precursors. These factors form a gene regulatory network which activates during LMPPs commitment towards B-lymphocytes development. We have assembled the functional interaction network, shown in Fig 1, by manually mining the information from the literature, as reported in S1 Text.

[Figure omitted. See PDF.]

Fig 1. The key TF interactions in the B-lymphopoiesis circuit.

The B-cell lineage is assumed to be determined by the antagonism between ZNF521 and the specification of the B-lymphoid lineage driven by E2A, EBF1 and PAX5. The interaction framework emerge out of literature. It is also supported by some assumptions from the phenomenological observation that can lead to cellular switching. In order to gain this particular behavior, considering well-known transcriptional motifs, we assumed a putative direct inhibition of ZNF521 through the regulator PAX5, (dashed line). Arrows represent gene activation, blunt arrows gene repression. TEBF1 and TZNF521, indicate environmental factors, which promote the transcriptional activity of EBF1 and ZNF521, respectively.

https://doi.org/10.1371/journal.pone.0132208.g001

From a dynamical system perspective, this biological network implements a switch-like function and, therefore, is expected to exhibit a bistable behavior, i.e., it can robustly operate in two distinct conditions. To enable the application of the analysis tools for bistability, we have complemented the qualitative description provided by the network in Fig 1 with kinetic rules associated to each interaction, thus devising the dynamical model Eqs (1a)–(1j). The model implicitly assumes that the activity levels of transcription factors and cytokine receptors are roughly proportional to their mRNA levels. In order to analyze their effect on early B-lymphopoiesis lineage commitment the model includes two environmental factors, TEBF1 and TZNF521, which promote the transcriptional activity of EBF1 and ZNF521, respectively. Modulating the effect of these exogenous inputs, we can induce the switch between the two stable operative conditions, analyzing the effect of EBF1 and ZNF521 on a continuous range of expression levels.

The dynamics of the system cannot be determined without first specifying the parameter values, fixing a point in the parameter space. The qualitative behavior of the commitment switch in the B lymphopoiesis lineage is adapted, firstly, by a set of biochemical parameters for the system in Eqs (1a)–(1j) that is reasonably consistent with the expression levels of microarray gene profiles of the CD19+ CD127+ and CD19+ CD127− cells with B lineage cell subsets isolated from pediatric marrow reported in [20, 21]. Since these earlier studies have been conducted using the HG-U133 chip set (HG-U133A and HG-U133B), analysis has been conducted over the 44,754 probe sets present in the two array types, which were then processed using the Affymetrix platform Expression Console. The application uses the Robust Multichip Analysis (RMA) in order to obtain the generation of experimental expression datasets. Further statistical analysis is performed considering Affymetrix Transcriptome Analysis Console (TAC) Software. Considering the different developmental lineages, we obtained a list of differentially expressed genes in genome-wide comparison between multipotent stage and B lymphoid LRPs (pro-B), as reported in Fig 2 and S1 Fig.

[Figure omitted. See PDF.]

Fig 2. Gene level differential expression analysis.

Every point in the scatter plots shows the expression of a gene in the two conditions CD34+/lin− vs. pro-B. HG-U133A Chip. Total number of genes evaluated: 22645, with 550 genes differentially expressed. CD34+/lin− vs. pro-B analysis shows that 167 genes are up-regulated and 383 genes are down-regulated. HG-U133B Chip. Total number of genes evaluated: 22283 with 1080 genes differentially expressed CD34+/lin− vs. pro-B analysis shows that 553 genes are up-regulated and 527 genes are down-regulated. According with the TAC default filter criteria, we setted a fold change (linear) cut-off of ±2 (all that points having a fold-change less than 2 are shown in gray) and p-value < 0,05. Red: up in CD34+/lin− vs pro-B; green: down in CD34+/lin− vs pro-B. Highlighted are shown the factors involved in the proposed model.

https://doi.org/10.1371/journal.pone.0132208.g002

Cell lineage determination in the state space

The state of our dynamical system, let us denote it by X, is defined by the expression levels of the genes included in the regulatory network in Fig 1, that is X = (xikaros, xgfi1, xpu.1, …, xcd19), where xikaros denotes the concentration of transcribed IKAROS mRNA. A point in the state space is a stable equilibrium if every trajectory of the system starting in a small neighborhood of that point remains close to it for an infinite time interval in the absence of exogenous inputs.

In order to implement a switch-like behavior, our system should exhibit two distinct equilibrium points: some genes will have high expression in one point and low in the other one, while other genes may settle on an intermediate expression level in both conditions. These two equilibrium points of our model represent the multipotent progenitors and committed B precursors conditions, and the model trajectories are attracted on one or the other depending on the initial condition [22]. In general, all the state trajectories are limited by gene regulatory interactions that represent effectively constraints in the dynamic evolution [23]. Indeed cellular development creates only a small subset of states among all the possible states because a cell state does not move entirely randomly. Obviously, individual genes do not alter their expression value independently because of predetermined regulatory interactions. No gene regulatory interactions entails free movement of states. In the proposed model, the admissible state trajectories are characterized by implicit constraints: for instance, if the transcription factor ZNF521 inhibits the expression of EBF1, then as ZNF521 increases its expression EBF1 is subjected to decrease.

A key circuit promotes irreversible bistable commitment switch between LMPPs and lymphoid LRPs

In correspondence to the lymphoid multipotent stage, the two transcription factors IKAROS and PU.1 initialize the lymphoid lineage commitment process, in particular activating the transcription of FLT3, a receptor specifically expressed on the surface of LMPPs [24, 25]. PU.1, in turn, is maintained at lower levels for the cell to persist in lymphoid line, avoiding the switching to the myeloid one. GFI1 (activated by IKAROS) has been suggested to interact in a direct manner with the promoter of SPI1 (the gene encoding PU.1 transcription), inhibiting its expression; hence, IKAROS indirectly lowers PU.1 level [26]. Our analysis has highlighted that the expression level of PU.1 cannot be simply nullified, but it must be finely regulated: indeed, its presence is necessary for the regulation of the expression of IL-7R, which, in turn, is required for the progression of the development from LMPPs to B-lymphoid lineage [27, 28].

The lymphoid LRPs commitment is guaranteed by the antagonism between ZNF521 and the B-lymphoid lineage specific genes, like TCF3 (coding for E2A trascription factor), EBF1 and PAX5, see Fig 3. While most of the interactions between these factors have been unraveled by means of biological experiments, evidence of regulation due to higher order multimeric bindings is lacking; therefore, in our model we admit only heterodimeric pairwise. If this modeling assumption was removed, one could rather easily (and arbitrarily) tailor high-order kinetics to obtain bistable dynamics. Our approach, instead, aims at the discovery of the key interactions underpinning the lineage commitment switch and are, therefore, the most promising candidates to be investigated via wet-lab experiments.

[Figure omitted. See PDF.]

Fig 3. Expression programs characterizing the commitment of LMPPs toward B-lymphoid LRPs.

Each state is guaranteed by a distinctive expression program of regulatory factors. In red (green), are depicted the factors up (down) regulated in each condition, LMPP and LRP, respectively. Values are taken in arbitrary units from the experimental gene profiling datasets, reported in [20].

https://doi.org/10.1371/journal.pone.0132208.g003

The model built on the basis of the currently known experimentally validated interactions is not able to reproduce the commitment switch behavior, i.e. it never achieves bistability, independently of the assigned parameter values. This suggests that part of the regulatory network still ought to be unveiled. Consequently, we have exploited the devised model to test additional hypothetical interactions, chosen among the most likely on the basis of the literature and of the knowledge of the domain expert.

The proposed in silico model, thus, provides an effective tool to rapidly sieve multiple hypotheses on the topology of the interaction network. From the examination of different network structures, only one interaction, namely the direct inhibition of the regulator PAX5 on ZNF521, has resulted to be crucial to reproduce the dynamics of a bistable lineage commitment switch. In fact, removing this inhibitory action of PAX5 on ZNF521 from the model, the system does no longer exhibit a switching characteristic, even in the presence of other additional interactions. For example, an alternative model featuring inhibition of ZNF521 through E2A cannot yield bistable behavior, see Model B in S1 Text and S6 Fig. Recently, mouse studies on the transcription factor Ikaros in early B-cell development have confirmed its role of regulator in lymphopoiesis, revealing also different activated or repressed target genes, among which the murine equivalent for ZNF521 [29]. Based on these experimental evidences, a further alternative model characterized by a direct repression on ZNF521 through IKAROS has been considered, Model C in S1 Text and S6 Fig. However, numerical simulations suggest that repression of IKAROS on ZNF521 is not able to induce bistable switching, see S6 Fig. Several other alternative interactions have been tested during the construction of the model, confirming that the PAX5/ZNF521 interaction is essential to the realization of the switching mechanism.

The mechanism linking EBF1 expression level to B-lymphoid LRPs commitment has the characteristics of an irreversible bistable switch, mediated by the feedback between ZNF521-EBF1, PAX5, and the multipotency marker FLT3. The system behavior changes are evaluated versus variations in control parameter regulating the transcriptional activity of EBF1 expression, TEBF1, by means of bifurcation analysis [30]. Externally modifying the value of TEBF1, it is possible to replicate experimental variation in the expression of EBF1, essential for the commitment of B-lymphocyte identity. Fig 4 shows a typical bistable characteristic, featuring the presence of two stable steady-state branches delimited by two bifurcation points. All the points included in the region of bistability admit either a low or high value in correspondence to a single value of the bifurcation parameter. Note also that the bifurcation point located on the negative semiaxis is inaccessible (expression levels cannot achieve negative values), therefore the two branches of the bifurcation diagram are disconnected, which renders the transition irreversible. In the bifurcation diagrams one branch corresponds to the multipotent progenitor condition and the other to the B-lymphoid committed progenitor condition. Once EBF1 level accumulates beyond the threshold, an irreversible transition occurs from multipotent progenitor (LMPP cell, stable state I) to B-lymphoid LRP (pro-B cell, stable state II).

[Figure omitted. See PDF.]

Fig 4. Irreversible and bistable switch in LMPPs commitment into B-lymphoid LRPs (pro-B cells).

Steady state concentrations of FLT3, ZNF521, PAX5 and CD19 as functions of an external signal TEBF1, which promotes B-cell development and represses ZNF521. All the points located on the negative semiaxis, the shaded area, are inaccessible (expression levels of some species assume negative values). LP, Limit Point.

https://doi.org/10.1371/journal.pone.0132208.g004

Analysis of the core feedback loops

In order to gain a deeper insight into the regulatory mechanisms that generate the lineage commitment, in the following we explore the dynamics of the different subnetworks that compose the core feedback loops. The expression and activity of EBF1 must be finely regulated in order to maintain the commitment of B-lymphoid development and balancing B lymphopoiesis [31, 32]. Among the inhibitors of EBF1, ZNF521 has been suggested to play a crucial role suppressing EBF1 activity [17, 18, 33]. More recent experiments have also shown that Ebf1 cooperates in a negative feedback loop to repress Zfp521 as differentiation proceeds [34]. Based on these experimental evidences, we examine the simplest feedback module, comprising only EBF1 and ZNF521, assuming heterodimeric bindings and positive autoregulations, S7 Fig (left panel). The associated mathematical model takes the form reported in Eq (2). Considering the nominal parameter values (S1 Table), the nullclines dZNF521/dt = 0 and dEBF1/dt = 0 intersect at a single steady state (SS) corresponding to the LMPP point, S7 Fig (upper right panel). The bifurcation diagrams in S7 Fig report the steady-state values of EBF1 and ZNF521 as function of the environmental factor TEBF1. These diagrams show the existence of a single inaccessible turning point (LP), thus system Eq (2) cannot provide a bistable switch. On the basis of this analysis, other molecular interactions have to be included in the model in order to to obtain bistable dynamics.

The commitment of LMPPs to LRPs is essentially sustained by the three transcription factors E2A, EBF1 and PAX5 [35]. E2A acts upstream of EBF1 by modulating its expression and hence, in concert with EBF1, activates the transcription of the PAX5 gene, leading to the reinforcement of the B-lymphoid commitment [36, 37]. Furthermore, EBF1 binds to and regulates in a direct manner the expression of PAX5 [38] which, in turn, stimulates the EBF1 expression through a positive feedback loop by binding to the proximal EBF1 promoter [39, 40]. Looking at the dynamical behavior of the EBF1/E2A/PAX5 positive feedback loop, see S8 Fig (left panel), it is possible to note that, although the cooperativity is increased by means of heterodimeric bindings between the TFs EBF1, E2A and PAX5, as described by Eq (3), the system cannot yield the classical S-shaped bifurcation plot that is characteristic of a bistable switch. Numerical continuation analysis, performed with respect to the environmental factor TEBF1, shows the existence of a single steady-state, suggesting that this feedback module fosters persistence in the committed state, S8 Fig (right panel).

Flt3-signaling cascade plays a very important role in B-cell development by priming early B-cell progenitors in order to proceed toward the B-cell differentiation, promoting the expression of IL-7 [41, 42]; indeed, mice lacking both FLT3- and IL-7R-derived signals fail to develop any B-cells [43]. Furthermore, the presence of a negative feedback on FLT3 transcription by PAX5 has been documented [44]. Therefore, another important issue is whether the negative feedback loop EBF1/PAX5/FLT3/IL-7R, illustrated in S9 Fig (left panel), may contribute to bistability. Bifurcation analysis performed with respect to TEBF1 values shows a single positive steady state for system Eq (4), relative to the fully committed state, in correspondence of which the expression of EBF1 and PAX5 are sustained at higher level and the FLT3 and IL-7R at lower ones, S9 Fig (right panel).

Expanding the model by assembling together the latter two feedback loops (see S10 Fig left panel) results in a reinforcement of the committed state, which does not provide the sought-after switch-like behavior, as derived by means of bifurcation analysis of model Eq (5), see S10 Fig (right panel). Analysis of the network obtained by combining the ZNF521-EBF1 loop with the EBF1/PAX5/FLT3/IL-7R one, see S11 Fig (left panel), does not appear to guarantee sufficient conditions for bistable behavior (the corresponding model is given in Eq (6): bifurcation analysis show the existence of positive steady-states relative to the committed lineage only, see S11 Fig (right panels). Also the combination of the feedback loop EBF1-ZNF521 with the EBF1/E2A/PAX5 one, see S12 Fig (left panel) and model Eq (7), does not yield bistability: bifurcation diagrams show that the system converges toward committed state.

Finally, we examine the dynamics generated by the hypothesized novel interaction, PAX5-ZNF521. Such interaction establishes a second feedback loop between EBF1 and ZNF521 mediated by PAX5, as illustrated in S13 Fig (left panel). This feedback structure, described by model Eq (8), is responsible for the irreversible commitment switch: bifurcation diagrams in S13 Fig (right panel) show typical curves exhibiting bistable irreversible behavior against variations in the transcriptional activation of EBF1.

Model assessment with respect to experimental observations

Similarly to other models of switch-like biological phenomena [11, 12], the proposed model is not aimed at providing quantitatively exact predictions of molecular concentrations, but rather to a system-level understanding of the mechanisms underlying the transitions between different operative conditions. Therefore, an assessment of the model plausibility can be done by means of a qualitative comparison of the in silico behavior with respect to the experimental bistable characteristics, rather than relying on some goodness-of-fit metrics. Thus, we compared the model steady-states with gene expression profiling datasets, focusing on the qualitative agreement between the model results and the gene profiling expressions (differentially gene expression results are reported in Fig 2 and S1 Fig). To confirm the phenomenological plausibility of the proposed model, we analyzed the results with respect to known molecular phenotypes, characterized by particular expression programs. We explicitly assessed the states of the attractors with gene expression profiles of differentiating multipotent progenitors toward B-lymphocytes. Human CD34+/lin− and precursor B-cell subsets specific expression values were taken from [20, 21]. Fig 5 reports a comparison between experimental expression data and the model equilibrium points, highlighting only those factors whose expression varies by a significant amount across the different cellular states, namely, EBF1, PAX5 and ZNF521, and the markers FLT3 and CD19, expressed on the cell surface in the multipotent and committed stage, respectively.

[Figure omitted. See PDF.]

Fig 5. Comparison between experimental expression profiling and computed attractors.

Each condition, for both experimental and model results, is normalized with respect to the given maximum value. It is possible to note a good agreement, in terms of bistable behavior, between the model prediction, the measured expression of the key factors involved in the commitment switch and the markers of the multipotent and the committed stage (FLT3 and CD19), respectively. Blue bars correspond to the LMPP stage, green bars to the LRP in the pro-B phase.

https://doi.org/10.1371/journal.pone.0132208.g005

For instance, in correspondence to the multipotent stage, the proposed model properly captures a higher expression of the LMPP phase promoters FLT3 and ZNF521; at the same time, the expression of the promoters of the lineage-restriction activity, EBF1, PAX5 and the surface B-precursor marker, CD19, is suppressed. In the B-lymphoid progenitors lineage, instead, the main fate B-cell development determinants EBF1, PAX5 and CD19 are strongly expressed in both experimental data and computational predictions, whereas FLT3 as ZNF521 are strongly repressed. These results suggest that the proposed commitment-switch model, derived by a combination of experimentally validated interactions and modeling hypotheses, is capable of reproducing the bistable behavior of the experimental system, in correspondence of the point in the parameter space defined by the solution proposed in S1 Table. This is not sufficient to determine the model plausibility: due to the significant inter-cellular variability of biomolecular circuits, a model can only be considered plausible if it is robust, i.e., the desired behavior is reproduced not only for a particular choice of the parameters, but over a non-trivial subset of the parameter-space [45–47].

Parametric sensitivity analysis showed that the bistable behavior of the proposed model is robust with respect to variations of the parameters involved in the dynamics of ZNF521 and EBF1. We have examined the range of parameter values in which bistable behavior is preserved, see S2 Table. In particular, the parameter b4, that expresses the strength of the inhibition of PAX5 on ZNF521 in the nonlinear ODE model Eqs (1a)–(1j), can assume a wide range of values, log2 fold change value: [-0.524; 2.944], computed as the log2 ratio between the minimum and maximum admissible value, respectively, for which the system is bistable and the parameter used in the model, S1 Table. Furthermore, we have been able to identify domains of variation of the parameters involved in the dynamics of IL-7R and E2A, which are characterized by bistable dynamics. This is in agreement with the fact that these factors are known to play a crucial role in B lymphopoiesis, as reported in the S1 Text. The parametric robustness of the proposed model supports its plausibility, suggesting that the qualitative response does not depend on the particular choice of the parameters, but is rather a structural feature emerging from the topology of the interaction network and the kinetics of these interactions.

Reprogramming strategies of B-lymphoid LRPs into multipotent progenitors

Directly increasing the EBF1 transcriptional activity, TEBF1, the EBF1 levels are sufficient to induce PAX5 expression, which, in turn, decreases ZNF521 and FLT3 levels, see Fig 4. Consequently, TCF3, EBF1 and PAX5 are activated, leading to lymphoid-LRP commitment (specified by expression of the surface marker CD19). The dynamics of these master regulators yield two stable regimes, corresponding to the two lymphoid differentiation stages. In principle, by modulating the external factors that regulate the expression of EBF1 and ZNF521, it is possible to revert the cell differentiation state from LRP to LMPP. To this aim, we need a strategy to make the system retrace back the classical S-shaped curves displayed in Fig 4. Since the bifurcations are driven by the master regulators EBF1, PAX5 and ZNF521, we have to act on these knobs to try to reprogram the cell fate.

In some cases it is sufficient to modify the expression of a single transcription factor within a gene regulatory network to reprogram the cell fate [48]. An important question is whether the mutually inhibitory feedback from EBF1 to ZNF521, and the feedback from PAX5 to FLT3 are essential to provide reinforcement of the commitment decision, and, more importantly, under which conditions the network can be modified to reverse the commitment. According to our model, the system can be reprogrammed from a state of lineage-restricted commitment by acting on the sole expression of ZNF521: this is evidenced by the bifurcation diagrams reported in Fig 6, showing that an increase in the expression of ZNF521 transforms the irreversible switch into a reversible one. Under unperturbed conditions, the irreversibility of the switch is guaranteed by the repressing action of PAX5 on ZNF521, which is directly influenced by the expression level of EBF1 but not by that of ZNF521. Therefore, to make the cell return to the preceding differentiation step, it is necessary to counteract this repression action by means of an exogenous transcriptional activation of ZNF521. As a consequence, the cell recovers the multipotent characteristics as self-renewal and differentiation into different lineages proceeding from a multipotent stage.

[Figure omitted. See PDF.]

Fig 6. EBF1-dependence of the B-cell development bifurcation diagram w.r.t the ZNF521 transcriptional activation.

(a) EBF1 versus external factor TEBF1 with no activation of ZNF521, TZNF521 = 0. Overexpression of EBF1 is sufficient to switch the system from the LMPP state (low EBF1) to the pro-B (high EBF1) state, as indicated by the arrows. (b) EBF1 versus external factor TEBF1 with transcriptional activation of ZNF521, TZNF521 = 0.12. A couple of limit points (LPs) defines a region of bistability for the proposed network. In this case, transcriptional activation of ZNF521 may convert the system from an irreversible to a reversible bistable switch. (c) EBF1 versus external factor TEBF1 with transcriptional activation of ZNF521, TZNF521 = 0.3. Consistently, higher values of transcriptional activation of ZNF521 entailing a sharp delay in B-cell development, showing a more ultrasensitive response to higher values of TEBF1.

https://doi.org/10.1371/journal.pone.0132208.g006

Discussion

Unveiling the mechanisms underpinning multipotent cell lineage commitment is one of the main goals of systems biology. The heterogeneity of multipotent and progenitor cells produces sub-populations that express slightly different properties and propensities to commit towards the various lineages.

In hematopoiesis, cell fate mechanisms at the level of multipotent progenitors commitment toward B-lymphoid progenitors are still only partially elucidated and, to our knowledge, no model has so far been developed to describe the regulation of this differentiation stage. We have described these mechanisms by means of a phenomenological mathematical model, to quantitatively characterize the transcriptional dynamics regulating the commitment of multipotent lymphoid progenitors to lymphoid LRPs. The proposed model focuses on the role of the transcription co-factor ZNF521, a stem cell-associated protein that has been implicated in the regulation of the homeostasis of the immature compartment in a variety of normal and neoplastic tissues [17, 18, 34, 49–57]. In particular, previous work from our group suggested that ZNF521 can antagonize B-cell development and contribute to maintain the multipotency of primitive lympho-myeloid progenitors, by counteracting their EBF1-driven B-lymphoid commitment [18]. In the present study, the different stages of lineage commitment and differentiation have been described mathematically in terms of stable attractors in the state-space of the devised dynamical model.

It is important to remark that the main goal of our work is to propose a phenomenological model useful for the comprehension of the regulatory mechanisms underlying some steps of B-lymphoid cell differentiation; for this reason, the model is not intended as a tool for quantitative prediction of gene expression but, rather, for describing the qualitative response of the system during the various differentiation stages and for pinpointing the key molecular mechanisms that enable the lineage commitment. Due to the lack of experimental evidences regarding the nature of the molecular reactions, only heterodimeric pairwise interactions have been considered in the design of the model. Although the restriction to this class of models limits the generality of our approach, it has been deemed necessary to a) limit the arbitrariness that is inherent in the process of deriving a model from a network of interactions, and to b) focus our study on the discovery of novel putative molecular mechanisms responsible for the lineage commitment switch. An alternative approach would have been to investigate whether the bistable behavior might be explained by the introduction of high-order cooperative effects, for example in the feedback loop EBF1-ZNF512. However, it is well known that by choosing a suitable complex kinetics it easy to induce a bistable behavior, even in very simple reaction networks (even with only two species). Therefore, the results obtained through this modeling approach would not have been useful as a discriminating factor for guiding future biological experiments.

The main objective of this work is to identify novel putative transcriptional interactions that mediate an irreversible switch-like commitment of multipotent progenitors toward the B-lymphoid lineage. ZNF521 has been shown to interact with multiple molecular partners, that are likely to mediate its regulatory activities [17, 34, 49–52, 55]. Among these, the inhibition of EBF1 appears to be of particular relevance in the maintenance of the homeostasis of the early hematopoietic cell compartment by opposing the activity of ZNF521[18].

A first contribution of the work has consisted in the design of a mathematical model that recapitulates the known regulatory interactions between the molecular factors involved in the lineage commitment of B-lymphoid cells, including ZNF521. Subsequently, this nominal model has been exploited to evaluate the plausibility of newly hypothesized interactions, not yet experimentally validated. These putative interactions have been translated into different extended models, which have been analyzed to establish their capability to exhibit the sought-after switch-like behavior. Model assessment has also been performed by comparison between the simulated expression profiles and experimental gene expression datasets, showing a good qualitative agreement between the simulated and the experimentally observed behavior.

Model-based quantitative analysis has shown that the currently known experimentally validated interactions are not sufficient to explain the irreversible commitment dynamics; starting from this finding. Among the various hypotheses that have been investigated, the most conceivable resulted to be the existence of an inhibitory effect of PAX5 on the transcription of ZNF521. PAX5 encodes a transcriptional regulator that has been shown to enable the B-cell-specific genetic program while repressing alternative ones [58]. Interestingly, Busslinger and colleagues have identified Pax5 binding within the Zfp521 locus in pro B cells [59], although they have not evidenced regulation of Zfp521 by Pax5. Direct repression of Zfp521, which is also highlighted in our model as the consequence of a possible direct interaction, has been reported in adipocyte progenitors by Kang et al. [34]. Considering that the expression of EBF1 precedes that of PAX5, it can be hypothesized that the initial ZNF521 repression by EBF1 is sufficient to lower the expression levels of the former, such that a further reduction of its expression in response to PAX5 may be undetectable by the gene profiling method used. Endowing our model with the PAX5 ⊣ ZNF521 interaction has resulted in the appearance of the commitment switch dynamics associated to the changes in the expression level of EBF1, PAX5 and ZNF521.

Sensitivity analysis of the model parameters has allowed us to assess the robustness of the commitment switch with respect to inter-cellular parametric variability, thus supporting the plausibility of the proposed mathematical description of the regulatory interaction network. Moreover, sensitivity analysis has supported the identification of the key interactions driving the lineage commitment mechanism: the parameters involved in the dynamics of EBF1 and ZNF521, that are a0, a2, …, a8, a10, …, a12, b0, b2, …, b4, exhibit a substantial range of variation in which switch like behavior is preserved; a quantitative evaluation of such ranges of variation is illustrated in S3 Fig and in S2 Table. The parameters f0, …, f3, j0, j1, which are involved in the dynamics of IL-7R and E2A, can also determine a loss of the switch-like behavior if their values fall outside certain intervals; S4 Fig and S2 Table report the detailed analysis performed on this parameters and the bounds of the corresponding bistability intervals. Such analysis has suggested that these factors are critical for transitions between the two solutions. Instead, the remaining parameters, involved in the dynamics of GFI1, PU.1, IKAROS (only its basal expression has been considered) FLT3, CD19, PAX5 show unimodal behavior for single variation, i.e. each parameter value is characterized by a single CD19 homeostatic level starting point, S5 Fig.

Under nominal conditions, the commitment mechanism has resulted to be irreversible, as expected according to the experimental observations (see Fig 4). By exploiting the devised model, we have been able to design a possible reprogramming strategy from the LRP to LMPP stage, which may provide a valuable indication for future experimental investigations, (see Fig 6). The analyses conducted on our model suggest that a) increased transcription of ZNF521 transforms the characteristic response of the commitment switch from an irreversible to a reversible form and b) the ZNF521-overxpressing cells may then be reverted to a multipotent state by repressing the expression of EBF1. It is interesting to notice, in this regard, that Zfp521 is among the 36 HSC-associated genes that have been initially used by Riddell et al. [60] to reprogram committed mouse hematopoietic cells to induced HSCs. Although in subsequent experiments a more restricted cocktail of factors that did not include Zfp521 has proven sufficient for this reprogramming, it will be of interest to test whether the combined overexpression of ZNF521 and silencing of EBF1 may indeed be effective in inducing the generation of multipotent progenitors from committed human B-lymphoid progenitors.

Dynamical models and the associated analysis methods provide effective tools for the comprehension of complex biological systems and for the in silico design of novel approaches to control their evolution. The regulatory mechanisms underlying multipotent cells differentiation represent a highly rewarding application field for such tools, as proven by previous analogous works, e.g. [12, 15, 16]. While the necessary steps of wet-lab experimental validation (such as the demonstration of the direct repression of ZNF521 by PAX5 or the consequences of enforced expression of ZNF521 in B-lymphoid progenitors that are currently underway) cannot be replaced by computational modeling and simulations, the presented results represent an advancement towards the comprehension of some key steps of B lymphopoiesis and will hopefully prove to be helpful in streamlining future experimental activity to better elucidate the hematopoietic regulatory mechanisms.

Materials and Methods

Determination of the transcriptional regulatory network

The first step in the generation of the proposed model consisted in the determination, based on the literature, of the key factors involved in the commitment of B-lymphoid lineage, and how these interact. The development of B lymphocytes from HSCs is organized as a multi-step programming process, turning on B-cell specific genes and silencing the expression of others. The commitment is initiated in correspondence of the LMPP stage, where the two transcription factors IKAROS and PU.1 start the early cellular specification [61, 62]. IKAROS and PU.1 can directly regulate the expression of FLT3, which is a marker of the LMPP stage [24, 25]. Our study also focused on the role of ZNF521 as an additional indicator of the multipotent developmental stage. At a later stage, the commitment of the lymphoid progenitor to B-lymphoid LRPs is essentially driven by three transcription factors: E2A, EBF1, and PAX5 [35]. This new stage is identified by a high expression level of CD19 on the cell surface. In the S1 Text, we provide a detailed analysis of the relevant literature, reporting the most up-to-date and significant experimental evidences concerning the LMPP to B-lymphoid LRP commitment mechanism. Furthermore, this biological information is arranged into a synthetic formal description, represented by the interaction network depicted in Fig 1.

Mathematical models

Full Model.

The gene regulatory network in Fig 1 has been translated into a mathematical model consisting of a set of nonlinear ordinary differential equations. The network dynamics have been modeled using simple expressions (e.g. quadratic and Michaelis-Menten-like mathematical terms), similarly to what is done with other phenomenological models of transcriptional networks, e.g. [12, 15]. Furthermore, for the sake of simplicity, we considered the following interpretation: the two cytokine signaling cascades on FLT3/Ras and IL-7R/Jak/Stat are described as direct interaction on the respective target. Hence, our model consists of only transcription factors and their interactions in a genetic control network, involving the ZNF521. We assume that the concentrations are in dimensionless units and the kinetic constants are in units of s−1, and the Michaelis-Menten constants are dimensionless. The corresponding dynamical equations gives the following(1a)(1b)(1c)(1d)(1e)(1f)(1g)(1h)(1i)(1j)where xIKAROS, xGFI1, xPU.1, xFKT3, xIL-7R, xE2A, xZNF521, xEBF1, xPAX5, and xCD19 stand for the corresponding expression levels, the parameters a0, b0, c0, d0, e0, f0, g0, h0, i0 and j0 represent the basal production rate of each protein and the μi with i = 1, …, 10, define the degradation rates. Each protein is transcribed by RNA polymerase when it is bound either by one or more transcription factors. For each transcription factor involved in the model, we assume the level of protein approximately proportional to the mRNA levels. In Eqs (1h) and (1g) we assumed by TEBF1 and TZNF521 the environmental factors activating the transcriptional activity of EBF1 and ZNF521, respectively, in order to analyze the effects on early B lymphopoiesis lineage commitment, promoting in turn the activation of one of these key regulators.

EBF1/ZNF521 feedback.

The dynamical equations corresponding to the simple feedback between EBF1 and ZNF521, reported in S7 Fig, assume the form:(2a)(2b)

EBF1/E2A/PAX5 feedback.

The dynamical equations corresponding to the EBF1/E2A/PAX5 positive feedback, reported in S8 Fig, are given by:(3a)(3b)(3c)

EBF1/PAX5/FLT3/IL-7R feedback.

The dynamical equations corresponding to the feedback involving EBF1, PAX5, FLT3 and IL-7R, reported in S9 Fig, are expressed as follows:(4a)(4b)(4c)(4d)

EBF1/E2A/PAX5/FLT3/IL-7R feedback.

The dynamical equations corresponding to the feedback involving EBF1, E2A, PAX5, FLT3 and IL-7R, reported in S10 Fig, are given by:(5a)(5b)(5c)(5d)(5e)

ZNF521/EBF1/PAX5/FLT3/IL-7R feedback.

The dynamical equations corresponding to the module constituted by ZNF521/EBF1/PAX5/FLT3/IL-7R, reported in S11 Fig, assume the form:(6a)(6b)(6c)(6d)(6e)

ZNF521/EBF1/E2A/PAX5 feedback.

The dynamical equations corresponding to the module constituted by ZNF521/EBF1/E2A/PAX5, reported in S12 Fig, assume the form:(7a)(7b)(7c)(7d)

EBF1/PAX5/ZNF521 feedback.

The dynamical equations corresponding to the EBF1/PAX5/ZNF521 feedback, reported in S13 Fig, assume the form:(8a)(8b)(8c)

Numerical simulations and bifurcation analysis have been performed by using Matlab and the Matcont continuation package [63], respectively. In order to get a biologically plausible model, we tuned the values of the model parameters such that the simulated response of the model matched the qualitative behavior of the experimental system, according to the experimental data reported in a set of microarray gene experiments. Given a certain network topology, along with a tentative model of its dynamics and admissible intervals of its parameters, the Bifurcation Discovery Tool (BDT) allows the computation of parameter set values that give rise to bifurcations [64]. The parameter space is explored by means of a genetic algorithm, seeking Hopf bifurcations, turning points and bistable switches. In particular, starting from an ODE-based reaction network, BDT allows an user to choose the parameters to be searched, the admissible parameter ranges, and the nature of the bifurcation to be sought. As a result, the tool will return the parameter values for which the sought-after behavior is observed. We have identified a set of parameter values that correspond to a bistable behavior of the network and to a system response that is in agreement with experimental measurements.

Microarray expression profiling

We considered gene expression given as supplementary data in [20, 21]. In the latter works, the authors treated the development of precursor B-cells from hemopoietic stem cells toward differentiation through a number of stages in the bone marrow before their migration to the periphery as naive mature B lymphocytes. (ArrayExpress Database, accession no. E-MEXP-384). B-cells progenitors were obtained from bone marrow samples of healthy children (age 3–16). Lineage-depleted (lin−) human cells enriched for expressing CD34 (CD34+ cells) were obtained from umbilical cord blood, since it is impossible to obtain enough CD34+/lin− cells from bone marrow samples. In general, lin− cells identify all the stem and progenitor cells for which mature cell lineage markers expression is undetectable. We have performed gene expression profiling using Affymetrix Expression Console, applying a robust multi-array normalization, and Transcriptome Analysis Console for differential expression analysis. Expression Console and TAC Softwares can be freely downloaded from the Affymetrix website.

Supporting Information

[Figure omitted. See PDF.]

S1 Text. Supporting Text.

https://doi.org/10.1371/journal.pone.0132208.s001

(PDF)

S1 Fig. Volcano plot showing gene level differential expression between the two conditions CD34+/lin− and pro-B vs. significance value.

All the points having a fold change (linear) greater than 2 and p-value < 0,05 indicate points-of-interest characterized by both large-magnitude fold-changes as well as high statistical significance (-log10 of p-value). All that points having a fold change (linear) less than 2 and p-value greater than 0,05 are shown in gray. Red: up in CD34+/lin− vs pro-B; green: down in CD34+/lin− vs pro-B. Highlighted are shown the factors involved in the proposed model confirming a well agreement between experimental and computed results.

https://doi.org/10.1371/journal.pone.0132208.s002

(EPS)

S2 Fig. Time profile concentrations of FLT3, ZNF521, EBF1, PAX5 and CD19 for the two cellular stages, indicating the final steady state values.

(a) When the ZNF521 is activated, the LMPPs specified genes are highly expressed arresting the development of B cells. Differently, (b) when no factor stimulates the expression of ZNF521, the system is driven towards the lymphoid lineage fate, resulting in a higher expression of the B-lymphoid LRP specified genes.

https://doi.org/10.1371/journal.pone.0132208.s003

(EPS)

S3 Fig. CD19 bifurcation analysis with respect to the parameters involved in EBF1 and ZNF521 transcriptional dynamics.

The parameters associated to EBF1 transcription, a0, a2, …, a8, a10, …, a12, and all the parameters involved in ZNF521 transcription, b0, b2, …, b4 exhibit an interval where bistable solutions can be found.

https://doi.org/10.1371/journal.pone.0132208.s004

(EPS)

S4 Fig. CD19 bifurcation analysis with respect to the parameters involved in IL-7R, PU.1, IKAROS and E2A transcriptional dynamics.

The parameters associated to IL-7R transcription, f0, …, f3 and to E2A transcription, j0, j1 exhibit an interval where bistable solutions can be found. Individual variations of the kinetic parameters describing the PU.1 dynamics, g0, …, g2, the GFI1 dynamics, h0, h1 and the basal expression of IKAROS, i0, instead, do not give rise to bistability.

https://doi.org/10.1371/journal.pone.0132208.s005

(EPS)

S5 Fig. CD19 bifurcation analysis with respect to the parameters involved in PAX5, CD19 and FLT3 transcriptional dynamics.

The parameters associated to PAX5 transcription, c0, …, c3, to CD19 transcription d0, d1 and all the parameters involved in FLT3 transcription, e1, …, e6 (with the only exception of its basal expression, e0), do not give rise to bistability.

https://doi.org/10.1371/journal.pone.0132208.s006

(EPS)

S6 Fig. Comparison of the proposed model with different hypothetical configurations.

Bifurcation study proves the crucial role of the inhibition of PAX5 on ZNF521 expression to admit the existence of commiment switching between LMPP and pro-B stages in B-lymphocyte development. Starting numerical continuation from the LMPP stage, let us impose the conditions related to the case of reversible commitment switch, for TZNF521 = 0.12. Model A. As already seen, ZNF521 inhibition through PAX5 generates a bistable commitment switch between LMPPs (lower values of CD19 expression) and LRP cells (greater values of CD19 expression). Model B. Removing ZNF521 inhibition through PAX5, the effects of inhibition of ZNF521 through the factor E2A have been considered. In this case, the system is characterized by a single CD19 expression level, preventing the occurrence of switching behavior. Model C. Removing ZNF521 inhibition through PAX5, the effects of ZNF521 repression through IKAROS have been considered. Also in this case, CD19 reaches a single expression level.

https://doi.org/10.1371/journal.pone.0132208.s007

(EPS)

S7 Fig. Dynamical analysis of the EBF1/ZNF521 feedback module.

The EBF1-ZNF521 mutual inhibition, with autoregulatory interactions; the environmental factors, which promote the transcriptional activity of EBF1 and ZNF521 are indicated as TEBF1 and TZNF521, respectively. The intersection among the nullclines dZNF521/dt = 0 and dEBF1/dt = 0, with parameters given by S1 Table, identifies a single stable point (SS) corresponding to the LMPP point (upper right panel). Bifurcation analysis performed promoting the transcriptional activity of EBF1 shows the existence of a single inaccessible turning point (LP), for negative values of TEBF1 (lower right panels). Solid line denotes stable steady-states, dashed line denotes unstable ones.

https://doi.org/10.1371/journal.pone.0132208.s008

(EPS)

S8 Fig. Dynamical analysis of the EBF1/E2A/PAX5 feedback module.

The positive feedback loop constituted by EBF1, E2A and PAX5, with autoregulatory loop and the environmental factor, TEBF1, acting on EBF1 transcription. Bifurcation analysis shows the expression levels of the corresponding factor as functions of TEBF1. It is possible to note that this feedback sustains the committed state, admitting only one stable positive point (right panels).

https://doi.org/10.1371/journal.pone.0132208.s009

(EPS)

S9 Fig. Dynamical analysis of the EBF1/PAX5/FLT3/IL-7R feedback module.

As suggested by experimental findings, EBF1, PAX5, FLT3 and IL-7R form a negative loop (left panel). Bifurcation analysis performed with respect to TEBF1 values shows that this feedback admits a single positive steady state, relative to the committed state, in correspondence of which the expression of EBF1 and PAX5 are sustained at higher level and the FLT3 and IL-7R at lower ones (right panels). Solid line denotes stable steady-states, dashed line denotes unstable ones.

https://doi.org/10.1371/journal.pone.0132208.s010

(EPS)

S10 Fig. Dynamical analysis of the EBF1/E2A/PAX5/FLT3/IL-7R feedback module.

The feedback loop comprising the positive loop involving EBF1, E2A and PAX5 and the negative one constituted by EBF1, PAX5, FLT3 and IL-7R (left panel). Bifurcation diagrams show the expression levels of the involved factors as functions of TEBF1. It is possible to observe a reinforcement of the committed state, which does not provide the sought-after switch-like behavior (right panels).

https://doi.org/10.1371/journal.pone.0132208.s011

(EPS)

S11 Fig. Dynamical analysis of the ZNF521/EBF1/PAX5/FLT3/IL-7R feedback module.

The module constituted by the mutual inhibition between ZNF521 and EBF1, with autoregulatory loops, and the negative feedback formed by EBF1/PAX5/FLT3/IL-7R (left panel), does not appear to be sufficient to bistability behavior. Bifurcation analysis performed versus TEBF1 show the existence of only positive steady-states relative to the committed lineage (right panels). Solid line denotes stable steady-states, dashed line denotes unstable ones.

https://doi.org/10.1371/journal.pone.0132208.s012

(EPS)

S12 Fig. Dynamical analysis of the ZNF521/EBF1/PAX5/E2A feedback module.

The module constituted by the mutual inhibition between ZNF521 and EBF1, with autoregulatory loops, and the positive feedback loop that sustain the committed state, EBF1/E2A/PAX5 (left panel), are unable to capture the required switch-like behavior. Observing the bifurcation diagrams (right panels) for consistent values of TEBF1 the system converges to the committed state. Solid line denotes stable steady-states, dashed line denotes unstable ones.

https://doi.org/10.1371/journal.pone.0132208.s013

(EPS)

S13 Fig. Dynamical analysis of the EBF1/PAX5/ZNF521 feedback module.

The EBF1/ZNF521 mutual inhibition, with autoregulatory loops, is enriched by a further interaction through the key factor PAX5. Observing the bifurcation diagrams (right panels) with respect to the TEBF1 values, this feedback induces irreversible bistable behavior. Solid line denotes stable steady-states, dashed line denotes unstable ones.

https://doi.org/10.1371/journal.pone.0132208.s014

(EPS)

S1 Table. Parameter values used in the dynamical equations.

Set of parameter values in correspondence of which the main model exhibits a classic bistable behavior.

https://doi.org/10.1371/journal.pone.0132208.s015

(PDF)

S2 Table. Parameter sensitivity coefficients for bistability.

The parameter sensitivity for bistability is defined as the log2 fold change of the ratio between the max (min) parameter value for which the system is bistable and the parameter used in the model, S1 Table.

https://doi.org/10.1371/journal.pone.0132208.s016

(PDF)

S3 Table. Parameter values used in the model B and C bifurcation studies.

Parameters values used as different configurations of the original transcriptional regulatory network, considering reasonable regulatory gene hypotheses.

https://doi.org/10.1371/journal.pone.0132208.s017

(PDF)

Acknowledgments

This work was partially supported by the Italian PON project “STAYWELL” (PON04a2_00003).

Author Contributions

Conceived and designed the experiments: LS CC GM FA. Performed the experiments: LS CC. Analyzed the data: LS CC GM. Wrote the paper: LS CC GM FA.

Citation: Salerno L, Cosentino C, Morrone G, Amato F (2015) Computational Modeling of a Transcriptional Switch Underlying B-Lymphocyte Lineage Commitment of Hematopoietic Multipotent Cells. PLoS ONE 10(7): e0132208. https://doi.org/10.1371/journal.pone.0132208

References

1. Zhu J, Emerson SG (2002) Hematopoietic cytokines, transcription factors and lineage commitment. Oncogene, 21(21):3295–313. pmid:12032771

2. Hu M, Krause D, Greaves M, Sharkis S, Dexter M, Heyworth C, et al. (1997) Multilineage gene expression precedes commitment in the hemopoietic system. Genes Dev, 11(6):774–785. pmid:9087431

3. Cross MA, Enver T (1997) The lineage commitment of haemopoietic progenitor cells. Curr Opin Genet Dev, 7(5):609–613. pmid:9388776

4. Swiers G, Patient R, Loose M (2006) Genetic regulatory networks programming hematopoietic stem cells and erythroid lineage specification. Dev Biol, 294(2):525–540. pmid:16626682

5. Loose M, Patient R (2006) Global genetic regulatory networks controlling hematopoietic cell fates. Curr Opin Hematol, 13(4):229–236. pmid:16755218

6. Loose M, Swiers G, Patient R (2007) Transcriptional networks regulating hematopoietic cell fate decisions. Curr Opin Hematol, 14(4):307–314. pmid:17534154

7. Graf T (2002) Differentiation plasticity of hematopoietic cells. Blood, 99(9):3089–3101. pmid:11964270

8. Laiosa CV, Stadtfeld M, Graf M (2006) Determinants of lymphoid-myeloid lineage diversification. Annu Rev Immunol, 24:705–738. pmid:16551264

9. Roeder I, Glauche I (2006) Towards an understanding of lineage specification in hematopoietic stem cells: a mathematical model for the interaction of transcription factors GATA-1 and PU.1. J Theor Biol, 241(4):852–65. pmid:16510158

10. Huang S, Guo YP, May G, Enver T (2007) Bifurcation dynamics in lineage-commitment in bipotent progenitor cells. Dev Biol, 305(2):695–713. pmid:17412320

11. Chickarmane V, Peterson C (2008) A computational model for understanding stem cell, trophectoderm and endoderm lineage determination. PLoS one, 3(10):e3478 pmid:18941526

12. Chickarmane V, Enver T, Peterson C, (2009) Computational modeling of the hematopoietic erythroid-myeloid switch reveals insights into cooperativity, priming, and irreversibility. PLoS Comput Biol, 5(1):e1000268. pmid:19165316

13. Laurent M, Kellershohn N (1999) Multistability: a major means of differentiation and evolution in biological systems. Trends Biochem Sci, 24(11):418–422. pmid:10542403

14. Arpin C, Déchanet J, Van Kooten C, Merville P, Grouard G, Brière F, et al. (1995) Generation of memory B-cells and plasma cells in vitro. Science, 268(5211):720–722. pmid:7537388

15. Zhang Q, Bhattacharya S, Kline DE, Crawford RB, Conolly RB, Thomas RS, et al. (2010) Stochastic modeling of B lymphocyte terminal differentiation and its suppression by dioxin. BMC Syst Biol, 4, 40. pmid:20359356

16. Martínez MR, Corradin A, Klein U, Álvarez MJ, Toffolo GM, di Camillo B, et al. (2012) Quantitative modeling of the terminal differentiation of B-cells and mechanisms of lymphomagenesis. Proc Natl Acad Sci USA, 109(7):2672–2677. pmid:22308355

17. Bond HM, Mesuraca M, Carbone E, Bonelli P, Agosti V, Amodio N, et al. (2004) Early hematopoietic zinc finger protein (EHZF), the human homolog to mouse Evi3, is highly expressed in primitive human hematopoieticv ells. Blood, 103(6):2062–2070. pmid:14630787

18. Mega T, Lupia M, Amodio N, Horton SJ, Mesuraca M, Pelaggi D, et al. (2011) Zinc finger protein 521 antagonizes early B-cell factor 1 and modulates the B lymphoid differentiation of primary hematopoietic progenitors. Cell Cycle, 10(13):2129–2139. pmid:21593590

19. Adams B, Dorfler P, Aguzzi A, Kozmik Z, Urbanek P, Maurer-Fogy I, et al. (1992) PAX-5 encodes the transcription factor BSAP and is expressed in B lymphocytes, the developing CNS, and adult testis. Genes Dev, 6(9):1589–1607. pmid:1516825

20. van Zelm MC, van der Burg M, de Ridder D, Barendregt BH, de Haas EFE, Reinders MJT, et al. (2005) Ig gene rearrangement steps are initiated in early human precursor B-cell subsets and correlate with specific transcription factor expression. J Immunol, 175(9):5912–5922. pmid:16237084

21. Nodland SE, Berkowska MA, Bajer AA, Shah N, de Ridder D, van Dongen JJM, et al. (2011) IL-7R expression and IL-7 signaling confer a distinct phenotype on developing human B-lineage cells. Blood, 118(8):2116–2127. pmid:21680796

22. Huang S, Kauffman S (2009) in Complex gene regulatory networks-from structure to biological observables: cell fate determination. Encyclopedia of Complexity and Systems Science, eds. Meyers RA (Springer), pp.1180–1293.

23. Huang S (2010) Cell Lineage Determination in State Space: A Systems View Brings Flexibility to Dogmatic Canonical Rules. PLoS Biol, 8(5): e1000380. pmid:20520792

24. Yoshida T, Ng SYM, Zuniga-Pflucker JC, Georgopoulos K (2006) Early hematopoietic lineage restrictions directed by Ikaros. Nat Immunol, 7(4):382–391. pmid:16518393

25. Carotta S, Dakic A, D’Amico A, Pang SHM, Greig KT, Nutt SL, et al. (2010) The transcription factor PU.1 controls dendritic cell development and FLT3 cytokine receptor expression in a dose-dependent manner. Immunity, 32(5):628–641. pmid:20510871

26. Spooner CJ, Cheng JX, Pujadas E, Laslo P, Singh H (2009) A recurrent network involving the transcription factors PU.1 and Gfi1 orchestrates innate and adaptive immune cell fates. Immunity, 31(4):576–586. pmid:19818654

27. De Koter RP, Lee HJ, Singh H (2002) PU.1 regulates expression of the interleukin-7 receptor in lymphoid progenitors. Immunity, 16(2):297–309.

28. Dakic A, Metcalf D, Di Rago L, Mifsud S, Li Wu, Nutt SL (2005) PU.1 regulates the commitment of adult hematopoietic progenitors and restricts granulopoiesis. J Exp Med, 201(9):1487–1502. pmid:15867096

29. Schwickert TA, Tagoh H, Gultekin S, Dakic A, Axelsson E, Minnich M, et al. (2014) Stage-specific control of early B cell development by the transcription factor Ikaros. Nat Immunol, 15:283–293. pmid:24509509

30. Strogatz SH (1994) Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering (Studies in Nonlinearity). Westview Press.

31. Lukin K, Fields S, Hartley J, Hagman J (2008) Early B-cell factor: Regulator of B lineage specification and commitment. Semin Immunol, 20(4):221–227. pmid:18722139

32. Nechanitzky R, Akbas D, Scherer S, Györy I, Hoyler T, Ramamoorthy S, et al. (2013) Transcription factor EBF1 is essential for the maintenance of B cell identity and prevention of alternative fates in committed cells. Nature Immunology, 14(8):867–875. pmid:23812095

33. Gerrits A (2011) in Exploiting Natural and Induced Genetic Variation to Study Hematopoiesis, PhD Thesis, University Library of Groningen.

34. Kang S, Akerblad P, Kiviranta R, Gupta RK, Kajimura S, Griffin MJ, et al. (2012) Regulation of early adipose commitment by Zfp521. PLoS Biol, 10(11):e1001433. pmid:23209378

35. Schebesta M, Heavey B, Busslinger M (2002) Transcriptional control of B-cell development Curr Opin Immunol, 14(2):216–223. pmid:11869895

36. Dias S, Månsson R, Gurbuxani S, Sigvardsson M, Kee BL (2008) E2A proteins promote development of lymphoid-primed multipotent progenitors. Immunity, 29(2):217–227. pmid:18674933

37. Kwon K, Hutter C, Sun Q, Bilic I, Cobaleda C, Malin S, et al. (2008) Instructive role of the transcription factor E2A in early B lymphopoiesis and germinal center B-cell development. Immunity, 28(6):751–762. pmid:18538592

38. Decker T, Pasca di Magliano M, McManus S, Sun Q, Bonifer C, Tagoh H, et al. (2009) Stepwise activation of enhancer and promoter regions of the B-cell commitment gene Pax5 in early lymphopoiesis. Immunity, 30(4):508–520. pmid:19345119

39. Roessler S, Györy I, Imhof S, Spivakov M, Williams RR, Busslinger M, et al. (2007) Distinct promoters mediate the regulation of Ebf1 gene expression by interleukin-7 and Pax5. Mol Cell Biol, 27(2):579–594. pmid:17101802

40. Fuxa M, Busslinger M (2007) Reporter gene insertions reveal a strictly B lymphoid-specific expression pattern of Pax5 in support of its B-cell identity function. J Immunol, 178(12):8222–8228. pmid:17600970

41. Miller JP, Izon D, DeMuth W, Gerstein R, Bhandoola A, Allman D (2002) The earliest step in B lineage differentiation from common lymphoid progenitors is critically dependent upon interleukin 7. J Exp Med, 196(5):705–711. pmid:12208884

42. Corfe SA, Paige CJ (2012) The many roles of IL-7 in B-cell development; mediator of survival proliferation and differentiation. Semin Immunol, 24(3):198–208. pmid:22421572

43. Sitnicka E, Brakebusch C, Martensson IL, Svensson M, Agace WW, Sigvardsson M, et al. (2003) Complementary signaling through flt3 and interleukin-7 receptor alpha is indispensable for fetal and adult B-cell genesis. J Exp Med, 198(10):1495–1506. pmid:14610045

44. Holmes ML, Carotta S, Corcoran LM, Nutt SL (2006) Repression of Flt3 by Pax5 is crucial for B-cell lineage commitment. Genes Dev, 20(8):933–938. pmid:16618805

45. Bates DG, Cosentino C (2011) Validation and invalidation of systems biology models using robustness analysis. IET Systems Biology, 5(4):229–244. pmid:21823754

46. Cosentino C, Bates DG (2011) Feedback Control in Systems Biology. CRC Press.

47. Kirk P, Thorne T, Stumpf MPH (2013) Model selection in systems and synthetic biology. Current opinion in Biotechnology, 24(4):767–774. pmid:23578462

48. Gurdon JB, Melton DA (2008) Nuclear reprogramming in cells. Science, 322(5909):1811–1815. pmid:19095934

49. Bond HM, Mesuraca M, Amodio N, Mega T, Agosti V, Fanello Delia, et al. (2008) Early hematopoietic zinc finger protein-zinc finger protein 521: a candidate regulator of diverse immature cells. Int J Biochem Cell Biol, 40(5):848–854. pmid:17543573

50. Matsubara E, Sakai I, Yamanouchi J, Fujiwara H, Yakushijin Y, Hato T, et al. (2009) The role of zinc finger protein 521/early hematopoietic zinc finger protein in erythroid cell differentiation. J Biol Chem, 284(6):3480–3487. pmid:19049973

51. Correa D, Hesse E, Seriwatanachai D, Kiviranta R, Saito H, Yamana K, et al. (2010) Zfp521 is a target gene and key effector of parathyroid hormone-related peptide signaling in growth plate chondrocytes. Dev Cell, 19(4):533–546. pmid:20951345

52. Hesse E, Saito H, Kiviranta R, Correa D, Yamana K, Neff L, et al. (2010) Zfp521 controls bone mass by HDAC3-dependent attenuation of Runx2 activity. J Cell Biol, 191(7):1271–1283. pmid:21173110

53. Yamasaki N, Miyazaki K, Nagamachi A, Koller R, Oda H, Miyazaki M, et al. (2010) Identification of Zfp521/ZNF521 as a cooperative gene for E2A-HLF to develop acute B-lineage leukemia. Oncogene, 29(13):1963–1975. pmid:20062079

54. Kamiya D, Banno S, Sasai N, Ohgushi M, Inomata H, Watanabe K, et al. (2011) Intrinsic transition of embryonic stem-cell differentiation into neural progenitors. Nature, 470(7335):503–509. pmid:21326203

55. Kiviranta R, Yamana K, Saito H, Ho DK, Laine J, Tarkkonen K, et al. (2013) Coordinated transcriptional regulation of bone homeostasis by Ebf1 and Zfp521 in both mesenchymal and hematopoietic lineages. J Exp Med, 210(5):969–985. pmid:23569325

56. Spina R, Filocamo G, Iaccino E, Scicchitano S, Lupia M, Chiarella E, et al. (2013) Critical role of zinc finger protein 521 in the control of growth, clonogenicity and tumorigenic potential of medulloblastoma cells. Oncotarget, 4(8):1280–1292. pmid:23907569

57. Mesuraca M, Galasso O, Leonardo G, Chiarella E, Scicchitano S, Vatrinet R, et al. (2014) Expression profiling and functional implications of a set of zinc finger proteins, ZNF423, ZNF470, ZNF521 and ZNF780B, in primary osteoarthritic articular chondrocytes. Mediators Inflamm, 2014:318793. pmid:24976683

58. Grosschedl R (2013) Establishment and Maintenance of B Cell Identity. Cold Spring Harb Symp Quant Biol, 157:549–564.

59. McManus S, Ebert A, Salvagiotto G, Medvedovic J, Sun Q, Tamir I, et al. (2011) The transcription factor Pax5 regulates its target genes by recruiting chromatin-modifying proteins in committed B cells. The EMBO Journal, 30(12):2388–2404. pmid:21552207

60. Riddell J, Gazit R, Garrison BS, Guo G, Saadatpour A, Mandal PK, et al. (2014) Reprogramming committed murine blood cells to induced hematopoietic stem cells with defined factors. Cell, 157:549–564. pmid:24766805

61. Kirstetter P, Thomas M, Dierich A, Kastner P, Chan S (2002) Ikaros is critical for B-cell differentiation and function. Eur J Immunol, 32(3):720–30. pmid:11870616

62. Nutt SL, Kee BL (2007) The transcriptional regulation of B-cell lineage commitment. Immunity, 26(6):715–725. pmid:17582344

63. Dhooge A, Govaerts W, Kuznetsov YA (2004) Numerical continuation of branch points of limit cycles in MATCONT. Lecture Notes in Computer Science, 3037:42–49.

64. Chickarmane V, Paladugu SR, Bergmann F, Sauro HM (2005) Bifurcation discovery tool. Bioinformatics, 21(18):3688–3690. pmid:16081475