 Proceedings
 Open Access
 Published:
Mathematical modeling of GATAswitching for regulating the differentiation of hematopoietic stem cell
BMC Systems Biologyvolume 8, Article number: S8 (2014)
Abstract
Background
Hematopoiesis is a highly orchestrated developmental process that comprises various developmental stages of the hematopoietic stem cells (HSCs). During development, the decision to leave the selfrenewing state and selection of a differentiation pathway is regulated by a number of transcription factors. Among them, genes GATA1 and PU.1 form a core negative feedback module to regulate the genetic switching between the cell fate choices of HSCs. Although extensive experimental studies have revealed the mechanisms to regulate the expression of these two genes, it is still unclear how this simple module regulates the genetic switching.
Methods
In this work we proposed a mathematical model to study the mechanisms of the GATAPU.1 gene network in the determination of HSC differentiation pathways. We incorporated the mechanisms of GATA switch into the module, and developed a mathematical model that comprises three genes GATA1, GATA2 and PU.1. In addition, a novel multipleobjective optimization method was designed to infer unknown parameters in the proposed model by realizing different experimental observations. A stochastic model was also designed to describe the critical function of noise, due to the small copy numbers of molecular species, in determining the differentiation pathways.
Results
The proposed deterministic model has successfully realized three stable steady states representing the priming and different progenitor cells as well as genetic switching between the genetic states under various experimental conditions. Using different values of GATA1 synthesis rate for the GATA1 protein availability in the chromatin sites during the time period of GATA switch, stochastic simulations for the first time have realized different proportions of cells leading to different developmental pathways under various experimental conditions.
Conclusions
Mathematical models provide testable predictions regarding the mechanisms and conditions for realizing different differentiation pathways of hematopoietic stem cells. This work represents the first attempt at using a discrete stochastic model to realize the decision of HSC differentiation pathways showing a multimodal distribution.
Background
Hematopoiesis is a highly orchestrated developmental process that comprises the proliferation, differentiation and maturation of a very small population of selfrenewing, pluripotent hematopoietic stem cells (HSC) for producing different types of blood cells, including erythrocyte, megakaryocyte, granulocyte, and macrophage [1, 2]. During development, the decision to leave the selfrenewing state and selection of a differentiation pathway are regulated by transcription factors (TFs) [3–6]. Intense experimental studies during the past two decades have suggested that tight regulation of HSC differentiation is controlled by the interaction of a number of genetic and epigenetic regulators of gene transcription, including the two TFs PU.1 and GATA1. Although the precise mechanisms to initiate the transcriptional cascade leading to different differentiated cells are not clear currently, experimental studies have established that both PU.1 and GATA1 'autoregulate' themselves, i.e. they stimulate their own production, as well as they are mutually antagonistic, i.e. they repress the production of each other [7–9]. In the erythrocyte/megakaryocite lineage high expression levels of gene GATA1 and low levels of PU.1 were detected [6, 10]; conversely, in the granulocyte/macrophage lineage higher expression levels of PU.1 and low levels of GATA1 were measured [5, 11]. However, the initial progenitor cells stay in the third state that has lowlevel activation of both genes PU.1 and GATA1. When a progenitor cell differentiates, it transitions from the initial indeterminate state into one of two differentiated states with high expression levels of either gene GATA1 or PU.1. Thus the full understanding of the PU.1 and GATA1 interaction is important for studying of the differentiation process of HSCs, which may be very useful in the clinical application of "differentiation therapy" for reestablishing the correct expression of PU.1 and/or GATA1 within immature leukemic cells [12].
To accurately describe the regulatory mechanisms controlling HSC differentiation, the important aspects that any mathematical model of the GATA1PU.1 network must include are an indeterminate state for the progenitor cells and two stable attractors of the dynamical system for the differentiated lineages. Since the first modeling attempt to study the regulation in the GATA1PU.1 network [13], a number of mathematical formalisms have been developed to realize the three stable steady states in HSCs. For example, Huang et al. used the Hill function with high cooperativity coefficients to qualitatively compare computer predictions with experimental evidence [14], giving support to the idea that lineage choice occurs as a two stage process, first priming and then differentiating. To remove the requirement of high cooperativity, Chickarmane et al. assumed that the autoregulation at both PU.1 and GATA1 occurs through the binding of monomers [15]. It was predicted that an additional mechanism should be involved in the repressive interaction to create a bistable switch, and therefore a third unspecified gene was introduced to create bistability. Alternatively, a mathematical model was proposed to include the dynamics of the inactive heterodimer GATA1PU.1 and the MichaelisMenten function was used to represent the low cooperativity [16]. We have designed a model by separating the strength of cooperativity for autoregulation and repression, and successfully realized a rich variety of systems behavior that have been found cross the existing models [17]. Recently Huang and collaborators proposed a general model that assumed no explicit interaction between the two genes; and realized a degenerated steady state via a new type of bifurcation [18]. Since the assumptions in these models are based on either unrealistic high cooperativity or unspecified gene, additional mechanisms are needed to accurately describe the differentiation decision of HSCs.
GATA2 is one of the six members of the hematopoietic GATA factor family and is most abundantly expressed in HSCs as well as in immature progenitors in hematopoietic lineages [10]. Previous studies demonstrated a critical role of GATA2 in the emergence and maintenance of HSCs [19]. In addition, strict regulation of genes GATA1 and GATA2 is critical for proper lineage commitment and development of erythroid cells. It was reported that GATA2 directly activated GATA1 expression in early erythroid progenitors, and then GATA1 accelerated its expression after its own expression has been initiated [20]. The gene expression process involving GATA1mediated displacement of GATA2 from chromatin is termed a GATA switch [21–23]. This "GATA factor switch" suggests a model whereby GATA2 and GATA1 sequentially bind the same cis elements with activating and repressive effects, respectively. During GATA factor switch, the GATA1 expression will increase due to the reciprocal decrease of GATA2, which leads progenitor cells to adopt an erythroid lineage. However, when GATA1 is absent or its expression is delayed, the reduction of GATA2 will increase the expression levels of PU.1 and leads progenitor cells to adopt a myeloid lineage [24, 25]. Therefore the dynamic expression patterns of GATA1 and GATA2 may influence the erythroidmyeloid cellfate selection by regulating the expression of gene PU.1 [22]. Although accumulating experimental evidence has suggested the important role of gene GATA2 in regulating the cellfate selection, only a simple Boolean network model has been proposed recently to include the regulatory function of gene GATA2 [26]. The kinetic dynamics of GATA2 and in particular the function of GATA switch has not been systematically studied so far.
Recent experimental studies have demonstrated that gene expression is a stochastic process. Key species of molecules such as DNA and mRNA may have small copy numbers, and the change of their copy numbers may cause significant variations of the system dynamics [27–30]. In particular, it has been shown that a variety of lineagerestricted genes in HSCs were expressed at low levels [31]. A recent study directly demonstrated that stochastic oscillation expression of lineageassociated genes could drive cellfate commitment [32]. Accumulating experimental evidence also suggested that stem cells are heterogeneous, with cells moving between two or more metastable states [33]. Recently a minimal model has been designed by combining cellextrinsic and cellintrinsic elements of regulation to understand how both instructive and stochastic events could inform cell commitment decision in hematopoiesis [34]. However, compared with the advances in developing various stochastic models to investigate the key functions of noise in genetic and cellular processes [35–37], the critical role of noise in determining the stem cell differentiation has not been well established. This work is aimed at developing the first stochastic model to explore the critical function of GATA switch and noise in determining the differentiation pathways of HSCs.
Methods
Mathematical model of GATAPU.1 gene network
In this work we propose a mathematical model for the GATAPU.1 regulatory network including genes GATA1, GATA2 and PU.1 (Figure 1A). It was assumed that each of these three genes activates their own expression through positive autoregulation. In addition, TF GATA2 activates the expression of gene GATA1 but inhibits the expression of gene PU.1, but the function of GATA2 to regulate the expressions of genes GATA1 and PU.1 is moderate [24]. Thus the expression levels of three genes GATA1, GATA2 and PU.1 may maintain at an intermediate state, which is compatible to the steady state of lineage priming. In addition, GATA1 inhibits the expression of GATA2 and PU.1, whereas PU.1 inhibits the expression of GATA1 and GATA2. Detailed assumptions of regulatory mechanisms and the list of chemical reactions are given in Additional file 1.
Based on the GATA switch model (Figure 1B), GATA2 localizes at the chromatin sites in early stage erythroblasts. When the expression levels of GATA1 increase as erythropoiesis progresses, GATA1 displaces GATA2 from chromatin sites and often (but not always) instigates a distant transcriptional output [21]. Although remaining in the cell, TF GATA2 in fact is unable to access the chromatin sites. To model this spatial regulatory mechanism, it was assumed that GATA2 proteins degrade during the process of GATA switch, which was realized by a large degradation rate constant ${k}_{2}^{*}$ of GATA2 in Eq. (1). Simultaneously, the concentration of GATA1 increases by using an additional synthesis rate of GATA1, which is proportional to the removal of GATA2.
Experiments also showed that, when GATA1 was not present, the shRNA knockdown of GATA2 increases the expression of gene PU.1 and reprograms the cells to become macrophages [24]. To realize genetic switching, the mechanism of the GATA switch and that of PU.1 knockdown were included in a single framework. Thus the knockdown of GATA2 was realized by a large degradation rate constant ${k}_{2}^{*}$ of GATA2 during a particular time period. In addition, a very small synthesis rate of GATA1 during that time period means the absence of GATA1 protein in the DNA promoter region. Furthermore, for simplicity, the basal expression rates of these three genes are assumed to be zero. We use the SheaAckers formalism, which is a widely used thermodynamic approach, to represent the gene expression based on the structure of transcription machinery [38]. All these assumptions led to the following model to realize the genetic switching of the GATAPU.1 regulatory network, given by
where x, y and z are the concentrations of TFs GATA1, GATA2 and PU.1, respectively, a_{1}, b_{1} and c_{1} represent the expression rates of genes GATA1, GATA2 and PU.1 autoregulated by itself, respectively, a_{2} is the expression rate of gene GATA1 regulated by TF GATA2, k_{1}, k_{2} and k_{3} are the degradation rates of TFs GATA1, GATA2 and PU.1, respectively, ${k}_{2}^{*}$ is the degradation rate constant to represent the displacement of GATA2 proteins,
during the time period [t_{1}, t_{2}] , and μ is a parameter to adjust the availability of GATA1 proteins in the chromatin sites. The GATA switch model is realized by using ${k}_{2}^{*}>0$, μ > 0 with a moderate value of parameter μ (e.g. μ = 1). In contrast, the knockdown of GATA2 is realized by using ${k}_{2}^{*}>0$, μ = 0. Another realization of the GATA2 knockdown is to use rate constants ${k}_{2}^{*}$ = 0, b_{1} = 0 for the expression of gene GATA2, namely to block the synthesis process of GATA2 but maintain the degradation process of GATA2 unchanged. Numerical results did not show any substantial difference between the computer simulations of these two types of realization (results not shown). Finally we note that the proposed model (Eq. 1) includes a recently designed model for the GATA1PU.1 module [18] as a special case if gene GATA2 is removed from the system.
Model parameters estimated from experimental data
There are 23 rate constants in the proposed mathematical model (Eq. 1). Part of the parameters was determined by the following experimental data.

(1)
The halflife of GATA1 protein is one hour [39]. In addition, the halflife of GATA2 is 30 min [40] that was confirmed by another observation in [33] therein. The halflife of PU.1 in Mel cell is ~2.4 h [41]. Thus the protein degradation rates constants are λ_{1} = 0.6931/h, λ_{2} = 1.3863/h, λ_{3} = 0.2888/h.

(2)
The disassociate rate of GATA1 binding to its DNA promoter is K_{ d } = 2.8 nM, which is more stable than the binding of GATA2 to its promoter site K_{ d } = 4.4 nM [42]. In addition, the disassociate rate of PU.1 binding to its DNA promoter is K_{ d } = 170 nM [43]. These rates were used to determine coefficients a_{4}, b_{4}, and c_{5}.

(3)
In addition, the heterodimer GATA1PU.1 has a 3fold increase of the binding rate constant over GATA1 to DNA [7]. It was assumed that a_{7} = 3a_{4}.

(4)
The disassociation rates of GATA1 and GATA2 binding to the DNA promoter sites are very close to each other [44]. It is assumed that the binding rates of GATA1 and GATA2 to the same DNA binding sites are the same, namely b_{6} = a_{7}, and c_{7} = c_{6}.
Multipleobjective optimization approach
To infer the unknown parameters in model (Eq. 1), we designed a novel multipleobjective optimization approach to estimate unknown model parameters (Figure 2). The criteria in this approach include the conditions of tristability, genetic switching and robustness property of the model. This approach includes the following major steps:
Step 1. Generate a set of model parameters in the GA. We first used the uniform distributed random variable U(0,1) generated in the genetic algorithm (GA) to determine the initial model rate constants. Note that the mutation manipulation in the GA was carried out on the random samples rather than model parameters. Since there are 7 restricted conditions for the existence of stable steady states (Eq. 6~8), other parameters were directly determined by the random samples a_{ i } = r_{ i }, where r_{ i } ~ U(0,1). Seven parameters were determined by the random samples together with the restricted conditions. For example, according to condition (Eq. 7), parameter b_{3} is determined by
where r_{10} ~ U (0,1) and ${x}_{1}^{*}$ is the steady state (Eq. 4). In addition, based on condition (Eq. 6), the four synthesis rates were also determined by a factor k,
where r_{ i } ~ U(0,1). We tested different values of k in order to realize genetic switching. When k is not large (1 < k < 100), we failed to find a set of parameters that could realize genetic switching. Thus in this work we used k = 1000 to search the unknown synthesis rate.
Step 2. Find the third steady state. The generated parameter set in Step 1 ensures the existence of the two stable steady states in which either GATA1 or PU.1 has high expression levels. To find the third stable steady state representing the primed progenitor state, we used MATLAB function fsolve.m to solving the nonlinear system for the steady state of the model (Eq. 1). If we can find the third stable steady state, we go to Step 3. Otherwise, we set the penalty function value to 4 and go to Step 6.
Step 3. Validate the existence of three steady states. Since the third steady state found in Step 2 was obtained by a numerical method, it may not exist due to the computational error. To examine the existence of the three steady states, we perturbed each steady state using samples of the uniformly distributed random variable and used the perturbed steady state as the initial condition to simulate system (Eq. 1). If the simulation converges to the steady state, it means the steady state exists, and then we go to Step 4. Otherwise, we set the penalty function value to 3 and go to Step 6.
Step 4. Validate the existence of genetic switching. We next examined whether the system can realize genetic switching using the mechanisms of GATA switch and GATA2 knockdown. The large degradation rate constant of GATA2 in model (Eq. 1) for these two mechanisms is
Simultaneously, the additional synthesis rate constant of GATA1 in model (Eq. 1) is
If the system realizes genetic switching through these two mechanisms, go to Step 5. Otherwise, set the penalty function value to 2 and go to Step 6.
Step 5. Robustness analysis. The inferred parameter set now satisfies the requirement for realizing three stable steady states and genetic switching. Next we used robustness property of the system as an additional criterion to choose the optimal model parameters. We tested the robustness property of mathematical model using the perturbed model parameters [45, 46]. Each parameter in the model was perturbed by a uniformly distributed random number
where U_{ jk } ~ U(0,1) and σ = 0.5 is the perturbation strength. For each set of model parameters, we obtained 1000 sets of perturbed parameters and then examined whether the mathematical model with the perturbed model parameters still maintained the three steady states. The model with a particular set of model parameters is more stable if the model maintains the three steady states with more sets of perturbed model parameters. To make an unbiased comparison, we used the same random samples U_{ jk } for different sets of parameters. The penalty function value at this step is the percentage of the parameter sets (from the 1000 sets of perturbed parameters) with which the model does not maintain the three steady states. Thus a smaller value of penalty function means better robustness property of the model.
Step 6. Return the value of penalty function to the genetic algorithm.
Note that there are four penalty functions in this proposed inference method, namely the existence of three steady states in Step 2 (denoted as event O_{1}), validation of these three steady states in Step 3 (denoted as O_{2}), validation of the existence of genetic switching in Step 4 (denoted as O_{3}), and robustness property of the model over 1000 perturbed sets of model parameters in Step 5 (denoted as O_{4}). Thus, if a numerical test in Steps 2~4 is not satisfied, the penalty score of that step should be larger than the maximal score of the next step. Since the maximal penalty score of Step 5 is set to unit one if none of the perturbed parameter set maintains genetic switching, the scores of the first three penalty functions are set to 4, 3 and 2, respectively. Using the notations of objective functions [47], the multiobjective function is represented by
Here the probability of objective O_{1} is defined by, for example,
Similar definitions are applied to the probabilities of objectives O_{2} and O_{3}.
Results
Steady states of the mathematical model
When ${k}_{2}^{*}$ = 0, the proposed mathematical model (Eq. 1) has up to four steady states. When k_{1}a_{4} ≠ 0 and k_{3}c_{5} ≠ 0, three steady states can be obtained analytically, given by
The existence conditions of these three steady states are given in Theorem 1.
Theorem 1 (1) The trivial steady state (Eq. 3) is unstable if any one of the following conditions is satisfied
(2) The steady state with high expression levels of gene GATA1 (Eq. 4) is stable if the following conditions are satisfied
(3) The steady state with high expression levels of gene PU.1 (Eq. 5) is stable if the following conditions are satisfied
The proof of this theorem is provided in the Additional file 1.
However, it is difficult to derive an analytical expression of the fourth potential steady state. For a given set of model parameters, we have to examine the existence of the fourth steady state numerically on a casebycase basis. The existence conditions of the stable steady states will be used as criteria to search the unknown model parameters.
Inference of model parameters
The proposed model (Eq. 1) has 20 unknown parameters by setting a_{3} = b_{2} = c_{2} = 1. We first estimated the values of 9 model parameters, namely (a_{4}, a_{7}, b_{4}, b_{6}, c_{5}, k_{1}, k_{2}, k_{3}) and (c_{6} = c_{7}) from the published experimental data discussed in the model section. Since there is not any published data for the temporal dynamics of gene expression, we used the designed multipleobjective optimization approach to estimate the remaining 11 unknown parameters (Figure 2). We used the genetic algorithm as an effective tool to search the optimal model parameters to realize genetic switching. A MATLAB toolbox developed by Chipperfield et al. [48] was used to infer the unknown model parameters. This toolbox used MATLAB functions to build a set of versatile routines for implementing a wide range of genetic algorithms. The initial estimate of rate constants can be changed by using different random seeds in the MATLAB toolbox, leading to different final estimates of the rate constants [49]. The genetic algorithm was run over 100 generations for each estimate, and we used a population of 1000 individuals in each generation. Our tests showed that a smaller population of individuals (i.e. 100 or 200) failed to produce an estimate of model parameters with which the model has tristability property. We implemented the genetic algorithm with different initial kinetic rates and obtained a number of estimated parameter sets. The parameter set having the best robustness property was selected as the final estimate.
Figure 3 gives simulations of the deterministic model (Eq. 1) using the degradation rate constant of GATA2 (Eq. 2) and different values of the GATA1 synthesis rate μ. The concentrations of GATA2 decreased substantially during t ∈ [500, 1500] when a large degradation rate ${k}_{2}^{*}$ = 6 was applied during that time period. Figures 3A, 3B and 3C give a simulation using the GATA switch module that was realized by μ = 1. When the concentration of GATA2 decreases during t ∈ [500, 1500] in Figure 3B, GATA1 reaches the high expression level quickly in Figure 3A. However, when GATA1 is absent (μ = 0) and GATA2 is knockdown, Figures 3D, 3E and 3F show a simulation leading to the high expression level of PU.1. The third scenario is that the expression of GATA1 maintains at a low level during the period of GATA2 decreases when μ = 0.32. In this case the expression levels of these three genes in Figures 3G, 3H and 3I maintain at an intermediate state that is dependent on the value of μ. Neither GATA1 nor PU.2 reached the steady state of high expression levels. When the expression level of GATA2 returned to the normal state at t = 1500, the system also returned to the primed steady state, which is an unsuccessful genetic switching.
When inferring the model parameters, we only tested the robustness property of the model with perturbation strength σ = 0.5 (see Step 5 in the Optimization Approach section). The next question is whether different values of the perturbation strength may lead to different selections of the estimated model parameters. To answer this question, we selected 10 sets of the estimated model parameters that have better robustness properties than the others, and examined the robustness property of the model using different perturbation strengths ranging from 0.1 to 1. Figure 4 shows that, the robustness property of the model with a smaller strength (σ = 0.3) is slightly different from that using the perturbation strength (σ = 0.5) because only a small numbers of perturbed parameters sets did not maintain the tristability property. However, when σ ≥ 0.4, the robustness property based on different values of perturbation strength σ is consistent with that using (σ = 0.5). Figure 4 also presents the averaged robustness properties based on 10 values of the perturbation strength (σ = 0.1, 0.2, ..., 1), which is well consistent with that using σ = 0.5. Thus it is suggested that the inferred model parameter is independent of the perturbation strength.
Bifurcation analysis
An important question in the estimation of model parameters is to what extent the model with the varied parameter still maintains the three stable steady states. To answer this question, we examined the tristability of the proposed model (Eq. 23) under the variations of one model parameter. For each rate constant, the perturbed values range from 0 to the 2fold of the estimated value. We first tested the influence of the four synthesis rates a_{1}, a_{2}, b_{1} and c_{1} since these parameters are important in regulating the expression of these three genes. Figure 5 shows that the tristability property is not sensitive to the change of synthesis rate a_{2}, which suggested that the expression of GATA1 is not sensitive to the regulation of GATA2. The result is consistent with the experimental observation showing that this regulation is relatively weaker than other genetic regulations in the system. Interestingly, the system maintains tristability if the synthesis rate of GATA1 (a_{1}) or PU.1 (c_{1}) increased, or the synthesis rate of GATA2 (b_{1}) decreased. However, the tristability property of the system disappears when these three rate constants vary along the opposite direction, in particular for the synthesis rates of GATA2 and PU.1. The discontinuity of the curves in Figure 5C suggested that the tristability property of the system disappears when the value of b_{1} is below a threshold value.
Stochastic dynamics
We have successfully realized the three steady states of the HSCs and genetic switching using the genetic regulations between GATA genes and gene PU.1. However, the proposed deterministic model failed to describe the heterogeneity in the differentiation decision of the HSCs. To tackle the challenge, we designed a stochastic model to realize the function of intrinsic noise arising from key species with small copy numbers, such as DNA and mRNA. Using our proposed stochastic modelling approach [50], we designed the following stochastic model based on the developed deterministic model (Eq. 1), given by
where X, Y, and Z are the copy numbers of protein GATA1, GATA2 and PU.1 respectively, $\text{}\tau $ is the stepsize, P(λ) is a Poisson random variable with mean λ. Protein concentrations in model Eq. (1) were transferred into the molecular copy numbers in this stochastic model. To reduce the computing time of stochastic simulations, the mechanisms of genetic switching were realized in the time period
Figure 6 gives three stochastic simulations using the same rate constants in Figure 4 and μ = 0.28. At time point t = 50, the expression level of GATA2 decreased quickly to the basal level due to the large degradation rate ${k}_{2}^{*}$. Although the expression levels of genes GATA1 and PU.1 increased quickly from time t = 50, the high expression level of one gene inhibited the expression of the other gene. Figures 6A, 6B and 6C show that the high expression levels of GATA1 inhibited the expression of gene PU.1 and GATA2, leading the cell into the erythrocyte lineage; whereas the expression levels of GATA1 failed to increase further in Figures 6D and its expression was inhibited by the growing expression levels of PU.1 in Figure 6F, leading the cell into the granulocyte lineage. In addition, it was possible that the expression levels of neither gene GATA1 nor PU.1 were high enough to inhibit the expression of the other gene. Figures 6G, 6H and 6I gave a simulation of unsuccessful switching in which the expression levels of these three genes maintained at an intermediate state that was determined by the value of μ. When the degradation rate of GATA2 returned to the normal value at t = 300, the expression levels of these three genes returned to the original priming state.
To obtain the statistical property of genetic switching, we generated 1000 stochastic simulations for each value of μ over the range of [0, 1.5]. When the value is small (μ < 0.11), which mimicked the experimental condition of gene GATA2 knockdown, Figure 7 shows that all simulations have high expression levels of gene PU.1. On the contrary, when the value of μ is large (μ ≥ 1.4), all simulations were switched to the state with high expression levels of gene GATA1, which represents the mechanisms of GATA switch that lead HSCs to the erythrocyte lineage. In addition, there is a range of values (0.12 <μ ≤ 0.22) with which the system failed to realize genetic switching and all stochastic simulations stay at the priming state. When the value of μ locates between these windows, stochastic simulations can realize two or three different steady states. For example, when the value of μ is 0.11 ≤ μ ≤ 0.12, the system state is in either the priming state or granulocyte lineage. An important prediction in Figure 7 is that, when the value of μ is (0.22 <μ ≤ 0.37), stochastic simulation of the system may lead to anyone of the three steady states, namely the erythrocyte, priming and granulocyte lineages, such as the case demonstrated in Figure 6.
Discussion
In this work we proposed a mathematical model to study the mechanisms of the GATAPU.1 gene network in the determination of HSC differentiation pathways. The novelty of this model is the inclusion of gene GATA2 and the GATA switch model based on the experimentally determined regulatory mechanisms. Our simulation results suggested that, based on the experimental determined regulatory mechanisms, the addition of the third gene, namely gene GATA2, is necessary and adequate to realize the three stable steady states of the HSCs. This result is consistent with the prediction that a connector gene X is required to realize the primed state [15]. Although the third gene negatively regulates the expression of PU.1 in these two models, the regulatory mechanisms of the third gene in these two models are not completely the same. For example, it was assumed that gene GATA1 positively regulate the expression of the unknown gene X [15], rather than the negative regulation of gene GATA2 by gene GATA1 in this work. Compared with the model in [15] in which additional and unknown external signals are required to maintain the tristability property, regulations between these three genes in our proposed model are adequate to maintain the tristability property of the system, which represents a successful approach in utilizing experimentally confirmed regulatory mechanisms to realize tristability property of the HSCs in regulating the erythroidmyeloid lineage decision.
The proposed model in this work for a network of three genes is a general framework that includes a recently published model for a network of two genes as a special case [18]. Unlike the model of two genes, which realized tristability using the bifurcation of the system by increasing the ratio of two model parameters, our model maintains the tristability property over a wide range of model parameter values. A related interesting question is the minimal motif to realize the stability property of a regulatory network with different numbers of steady states. It has been demonstrated that the twogene module with selfactivation and mutual repression can realize stability property with two steady states [50–53]. Although attempts have been made to realize tristability property of this twogene module using various assumptions of regulatory mechanisms [13, 14, 17], our research demonstrated that a regulatory module with three genes could maintain stability property with three steady states without the assumption of the autoregulation via high order multimers. A similar result is that a threecomponent motif with four links can realize tristability [54]. However, the GATAPU.1 module has six links together with more complex regulatory mechanisms. It would be interesting to analyze theoretically the stability property of a regulatory module with the maximal number of steady states under various regulatory mechanisms.
Stochastic simulations in this work predicted that the synthesis rate of GATA1 during the decreasing process of GATA2 determines the probability of erythroidmyeloid lineage decision. Since this synthesis rate represents the availability of GATA1 in the genomic regulatory regions during the GATA switch, there are two potential resources of noise to explain the variations of GATA1 proteins in the genomic regions. First, recent experimental research investigating cellular processes at single cells has revealed convincing evidence showing large heterogeneity in protein abundance and dynamics among genetically identical cells [55]. The variations of protein concentration in different cells may lead to different rates for GATA1 to enter the regulatory regions during the process of GATA switch. The heterogeneity in protein abundance may be one of the reasons to explain the differentiation preference of HSCs for the erythroid lineage or the myeloid lineage. In addition, the intrinsic noise due to the small copy numbers of molecular species in the GATAPU.1 module may further contribute the heterogeneity to the HSC lineage differentiation decision even for cells having the same lineage decision preference. These various resources of noise in gene expression may be part of the transcriptomewide noise proposed by Huang and colleagues [32]. Therefore more comprehensive stochastic models are needed to explain the functions of other types of noise in the decision of lineage selection.
Conclusions
In summary, we proposed a mathematical model to study the mechanisms of the GATAPU.1 gene network in the determination of HSC differentiation pathways. In addition, a multipleobjective optimization approach was developed to infer model parameters in order to realize the three stable steady states representing the three different types of blood cells and genetic switching. A stochastic model was also designed to describe the function of noise in determining the differentiation pathways. Stochastic simulations successfully realized different proportions of cells leading to different developmental pathways under the same experimental conditions, and provided testable predictions regarding the conditions and mechanisms to realize different differentiation pathways. This work represents the first attempt at using a discrete stochastic model to realize the decision of HSC differentiation pathways with a multimodal distribution.
References
 1.
Wang LD, Wagers AJ: Dynamic niches in the origination and differentiation of haematopoietic stem cells. Nat Rev Mol Cell Biol. 2011, 12 (10): 643655. 10.1038/nrm3184.
 2.
Cedar H, Bergman Y: Epigenetics of haematopoietic cell development. Nat Rev Immunol. 2011, 11 (7): 478488. 10.1038/nri2991.
 3.
Mercer EM, Lin YC, Murre C: Factors and networks that underpin early hematopoiesis. Semin Immunol. 2011, 23 (5): 317325. 10.1016/j.smim.2011.08.004.
 4.
Dore LC, Crispino JD: Transcription factor networks in erythroid cell and megakaryocyte development. Blood. 2011, 118 (2): 231239. 10.1182/blood201104285981.
 5.
Friedman AD: Transcriptional control of granulocyte and monocyte development. Oncogene. 2007, 26 (47): 68166828. 10.1038/sj.onc.1210764.
 6.
Goldfarb AN: Transcriptional control of megakaryocyte development. Oncogene. 2007, 26 (47): 67956802. 10.1038/sj.onc.1210762.
 7.
Liew CW, Rand KD, Simpson RJ, Yung WW, Mansfield RE, Crossley M, ProetoriusIbba M, Nerlov C, Poulsen FM, Mackay JP: Molecular analysis of the interaction between the hematopoietic master transcription factors GATA1 and PU.1. J Biol Chem. 2006, 281 (38): 2829628306. 10.1074/jbc.M602830200.
 8.
Nishimura S, Takahashi S, Kuroha T, Suwabe N, Nagasawa T, Trainor C, Yamamoto M: A GATA box in the GATA1 gene hematopoietic enhancer is a critical element in the network of GATA factors and sites that regulate this gene. Mol Cell Biol. 2000, 20 (2): 713723. 10.1128/MCB.20.2.713723.2000.
 9.
Okuno Y, Huang G, Rosenbauer F, Evans EK, Radomska HS, Iwasaki H, Akashi K, MoreauGachelin F, Li Y: Potential autoregulation of transcription factor PU.1 by an upstream regulatory element. Mol Cell Biol. 2005, 25 (7): 28322845. 10.1128/MCB.25.7.28322845.2005.
 10.
Akashi K, Traver D, Miyamoto T, Weissman IL: A clonogenic common myeloid progenitor that gives rise to all myeloid lineages. Nature. 2000, 404 (6774): 193197. 10.1038/35004599.
 11.
Gupta P, Gurudutta GU, Saluja D, Tripathi RP: PU.1 and partners: regulation of haematopoietic stem cell fate in normal and malignant haematopoiesis. J Cell Mol Med. 2009, 13 (1112): 43494363. 10.1111/j.15824934.2009.00757.x.
 12.
Burda P, Laslo P, Stopka T: The role of PU.1 and GATA1 transcription factors during normal and leukemogenic hematopoiesis. Leukemia. 2010, 24 (7): 12491257. 10.1038/leu.2010.104.
 13.
Roeder I, Glauche I: Towards an understanding of lineage specification in hematopoietic stem cells: a mathematical model for the interaction of transcription factors GATA1 and PU.1. J Theor Biol. 2006, 241 (4): 852865. 10.1016/j.jtbi.2006.01.021.
 14.
Huang S, Guo YP, May G, Enver T: Bifurcation dynamics in lineagecommitment in bipotent progenitor cells. Dev Biol. 2007, 305 (2): 695713. 10.1016/j.ydbio.2007.02.036.
 15.
Chickarmane V, Enver T, Peterson C: Computational modeling of the hematopoietic erythroidmyeloid switch reveals insights into cooperativity, priming, and irreversibility. PLoS Comput Biol. 2009, 5 (1): e100026810.1371/journal.pcbi.1000268.
 16.
Bokes P, King JR, Loose M: A bistable genetic switch which does not require high cooperativity at the promoter: a twotimescale model for the PU.1GATA1 interaction. Math Med Biol. 2009, 26 (2): 117132. 10.1093/imammb/dqn026.
 17.
Duff C, SmithMiles K, Lopes L, Tian T: Mathematical modelling of stem cell differentiation: the PU.1GATA1 interaction. J Math Biol. 2012, 64 (3): 449468. 10.1007/s0028501104193.
 18.
Andrecut M, Halley JD, Winkler DA, Huang S: A general model for binary cell fate decision gene circuits with degeneracy: indeterminacy and switch behavior in the absence of cooperativity. PLoS One. 2011, 6 (5): e1935810.1371/journal.pone.0019358.
 19.
Ling KW, Ottersbach K, van Hamburg JP, Oziemlak A, Tsai FY, Orkin SH, Ploemacher R, Hendriks RW, Dzierzak E: GATA2 plays two functionally distinct roles during the ontogeny of hematopoietic stem cells. J Exp Med. 2004, 200 (7): 871882. 10.1084/jem.20031556.
 20.
Grass JA, Boyer ME, Pal S, Wu J, Weiss MJ, Bresnick EH: GATA1dependent transcriptional repression of GATA2 via disruption of positive autoregulation and domainwide chromatin remodeling. Proc Natl Acad Sci USA. 2003, 100 (15): 88118816. 10.1073/pnas.1432147100.
 21.
Bresnick EH, Lee HY, Fujiwara T, Johnson KD, Keles S: GATA switches as developmental drivers. J Biol Chem. 2010, 285 (41): 3108731093. 10.1074/jbc.R110.159079.
 22.
Kaneko H, Shimizu R, Yamamoto M: GATA factor switching during erythroid differentiation. Curr Opin Hematol. 2010, 17 (3): 163168.
 23.
Snow JW, Trowbridge JJ, Johnson KD, Fujiwara T, Emambokus NE, Grass JA, Orkin SH, Bresnick EH: Contextdependent function of "GATA switch" sites in vivo. Blood. 2011, 117 (18): 47694772. 10.1182/blood201010313031.
 24.
Chou ST, Khandros E, Bailey LC, Nichols KE, Vakoc CR, Yao Y, Huang Z, Crispino JD, Hardison RC, Blobel GA: Graded repression of PU.1/Sfpi1 gene transcription by GATA factors regulates hematopoietic cell fate. Blood. 2009, 114 (5): 983994. 10.1182/blood200903207944.
 25.
Huang Z, Dore LC, Li Z, Orkin SH, Feng G, Lin S, Crispino JD: GATA2 reinforces megakaryocyte development in the absence of GATA1. Mol Cell Biol. 2009, 29 (18): 51685180. 10.1128/MCB.0048209.
 26.
Krumsiek J, Marr C, Schroeder T, Theis FJ: Hierarchical differentiation of myeloid progenitors is encoded in the transcription factor network. PLoS One. 2011, 6 (8): e2264910.1371/journal.pone.0022649.
 27.
Raser JM, O'Shea EK: Noise in gene expression: origins, consequences, and control. Science. 2005, 309 (5743): 20102013. 10.1126/science.1105891.
 28.
Raj A, van Oudenaarden A: Nature, nurture, or chance: stochastic gene expression and its consequences. Cell. 2008, 135 (2): 216226. 10.1016/j.cell.2008.09.050.
 29.
Balazsi G, van Oudenaarden A, Collins JJ: Cellular decision making and biological noise: from microbes to mammals. Cell. 2011, 144 (6): 910925. 10.1016/j.cell.2011.01.030.
 30.
Hume DA: Probability in transcriptional regulation and its implications for leukocyte differentiation and inducible gene expression. Blood. 2000, 96 (7): 23232328.
 31.
Miyamoto T, Iwasaki H, Reizis B, Ye M, Graf T, Weissman IL, Akashi K: Myeloid or lymphoid promiscuity as a critical step in hematopoietic lineage commitment. Dev Cell. 2002, 3 (1): 137147. 10.1016/S15345807(02)002010.
 32.
Chang HH, Hemberg M, Barahona M, Ingber DE, Huang S: Transcriptomewide noise controls lineage choice in mammalian progenitor cells. Nature. 2008, 453 (7194): 544547. 10.1038/nature06965.
 33.
Graf T, Stadtfeld M: Heterogeneity of embryonic and adult stem cells. Cell Stem Cell. 2008, 3 (5): 480483. 10.1016/j.stem.2008.10.007.
 34.
Palani S, Sarkar CA: Integrating extrinsic and intrinsic cues into a minimal model of lineage commitment for hematopoietic progenitors. PLoS Comput Biol. 2009, 5 (9): e100051810.1371/journal.pcbi.1000518.
 35.
Wilkinson DJ: Stochastic modelling for quantitative description of heterogeneous biological systems. Nat Rev Genet. 2009, 10 (2): 122133. 10.1038/nrg2509.
 36.
Kaern M, Elston TC, Blake WJ, Collins JJ: Stochasticity in gene expression: From theories to phenotypes. Nature Reviews Genetics. 2005, 6 (6): 451464. 10.1038/nrg1615.
 37.
Raj A, van Oudenaarden A: Nature, Nurture, or Chance: Stochastic Gene Expression and Its Consequences. Cell. 2008, 135 (2): 216226. 10.1016/j.cell.2008.09.050.
 38.
Shea MA, Ackers GK: The OR control system of bacteriophage lambda. A physicalchemical model for gene regulation. J Mol Biol. 1985, 181 (2): 211230. 10.1016/00222836(85)900865.
 39.
HernandezHernandez A, Ray P, Litos G, Ciro M, Ottolenghi S, Beug H, Boyes J: Acetylation and MAPK phosphorylation cooperate to regulate the degradation of active GATA1. Embo J. 2006, 25 (14): 32643274. 10.1038/sj.emboj.7601228.
 40.
Ferreira R, Wai A, Shimizu R, Gillemans N, Rottier R, von Lindern M, Ohneda K, Grosveld F, Yamamoto M, Philipsen S: Dynamic regulation of Gata factor levels is more important than their identity. Blood. 2007, 109 (12): 54815490. 10.1182/blood200611060491.
 41.
Papetti M, Skoultchi AI: Reprogramming leukemia cells to terminal differentiation and growth arrest by RNA interference of PU.1. Molecular cancer research: MCR. 2007, 5 (10): 10531062. 10.1158/15417786.MCR070145.
 42.
Ghirlando R, Trainor CD: Determinants of GATA1 binding to DNA: the role of nonfinger residues. J Biol Chem. 2003, 278 (46): 4562045628. 10.1074/jbc.M306410200.
 43.
Pio F, AssaMunt N, Yguerabide J, Maki RA: Mutants of ETS domain PU.1 and GGAA/T recognition: free energies and kinetics. Protein science: a publication of the Protein Society. 1999, 8 (10): 20982109. 10.1110/ps.8.10.2098.
 44.
Ko LJ, Engel JD: DNAbinding specificities of the GATA transcription factor family. Mol Cell Biol. 1993, 13 (7): 40114022.
 45.
Kitano H: Biological robustness. Nat Rev Genet. 2004, 5 (11): 826837. 10.1038/nrg1471.
 46.
Kitano H: Towards a theory of biological robustness. Mol Syst Biol. 2007, 3: 137
 47.
Wang Y, Wu QF, Chen C, Wu LY, Yan XZ, Yu SG, Zhang XS, Liang FR: Revealing metabolite biomarkers for acupuncture treatment by linear programming based feature selection. BMC Syst Biol. 2012, 6 (Suppl 1): S1510.1186/175205096S1S15.
 48.
Chipperfield AJ, Fleming PJ, Fonseca CM: Genetic Algorithm Tools for Control Systems Engineering. Proceedings of Adaptive Computing in Engineering Design and Control. 1994, Plymouth Engineering Design Centre, 128133.
 49.
Tian T, Xu S, Gao J, Burrage K: Simulated maximum likelihood method for estimating kinetic rates in gene expression. Bioinformatics. 2007, 23 (1): 8491. 10.1093/bioinformatics/btl552.
 50.
Tian T, Burrage K: Stochastic models for regulatory networks of the genetic toggle switch. Proc Natl Acad Sci USA. 2006, 103 (22): 83728377. 10.1073/pnas.0507818103.
 51.
Tian T, Burrage K: Bistability and switching in the lysis/lysogeny genetic regulatory network of bacteriophage lambda. J Theor Biol. 2004, 227 (2): 229237. 10.1016/j.jtbi.2003.11.003.
 52.
Gardner TS, Cantor CR, Collins JJ: Construction of a genetic toggle switch in Escherichia coli. Nature. 2000, 403 (6767): 339342. 10.1038/35002131.
 53.
Kobayashi H, Kaern M, Araki M, Chung K, Gardner TS, Cantor CR, Collins JJ: Programmable cells: interfacing natural and engineered gene networks. Proc Natl Acad Sci USA. 2004, 101 (22): 84148419. 10.1073/pnas.0402940101.
 54.
Tyson JJ, Novak B: Functional motifs in biochemical reaction networks. Annu Rev Phys Chem. 2010, 61: 219240. 10.1146/annurev.physchem.012809.103457.
 55.
Yuan TL, Wulf G, Burga L, Cantley LC: Celltocell variability in PI3K protein level regulates PI3KAKT pathway activity in cell populations. Curr Biol. 2011, 21 (3): 173183. 10.1016/j.cub.2010.12.047.
Acknowledgements
This work is supported by grants from the Australian Research Council (ARC) (DP1094181 and DP120104460). T.T. is also in receipt of an ARC Future Fellowship (FT100100748).
Declarations
The publication costs for this article were funded by the corresponding author.
This article has been published as part of BMC Systems Biology Volume 8 Supplement 1, 2014: Selected articles from the Twelfth Asia Pacific Bioinformatics Conference (APBC 2014): Systems Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcsystbiol/supplements/8/S1.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
T.T conceived the research. T.T. and K.SM. conducted research, analyzed data and wrote the manuscript. All authors read and approved the final manuscript.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Published
DOI
Keywords
 Stable Steady State
 Robustness Property
 Genetic Switching
 Perturbation Strength
 Chromatin Site