Original Research ARTICLE
Asymmetry in erythroid-myeloid differentiation switch and the role of timing in a binary cell-fate decision
- 1Nonlinear Analysis and Applied Mathematics Research Group (NAAM), Department of Mathematics, King Abdulaziz University, Jeddah, Saudi Arabia
- 2Department of Mathematics and Institute for Women’s Health, University College London, London, UK
GATA1-PU.1 genetic switch is a paradigmatic genetic switch that governs the differentiation of progenitor cells into two different fates, erythroid and myeloid fates. In terms of dynamical model representation of these fates or lineages corresponds to stable attractor and choosing between the attractors. Small asymmetries and stochasticity intrinsically present in all genetic switches lead to the effect of delayed bifurcation which will change the differentiation result according to the timing of the process and affect the proportion of erythroid versus myeloid cells. We consider the differentiation bifurcation scenario in which there is a symmetry-breaking in the bifurcation diagrams as a result of asymmetry in external signaling. We show that the decision between two alternative cell fates in this structurally symmetric decision circuit can be biased depending on the speed at which the system is forced to go through the decision point. The parameter sweeping speed can also reduce the effect of asymmetry and produce symmetric choice between attractors, or convert the favorable attractor. This conversion may have important contributions to the immune system when the bias is in favor of the attractor which gives rise to non-immune cells.
The importance of studying the immune system has attracted mathematicians and biologists to discover more of its features in recent years. One of the mechanisms is to study the genetic networks that control the lineage commitment of hematopoietic stem cells, which produce the full range of blood cells, including the immune cells (1). Many mathematical models have been used to study the differentiation of progenitor cell into erythroid and myeloid lineages based on the expression of lineage-specific transcription factors GATA1 and PU.1, respectively (2, 3). An important question arises in these models about the causes of bifurcation and symmetry-breaking and whether they occur in response to intrinsic cues or extrinsic signals. In fact, the integration of both intrinsic and extrinsic factors has received an extensive attention to elucidate the roles of external signals in cell-fate decision processes, and most importantly its relationship to the production of immune cells (3–6). Another important and interesting factor that can affect the decision of the cell is the speed of external signals or the speed of crossing the critical region (7–9). Remarkably, varying control parameter with time has been studied in many other systems. Ashwin et al. (10) have investigated how the rate of change of a parameter (or input) imposes significant changes in the climate system. It is found that rapid change may force the system to move away from a branch of attractors. This dependence on the rate was referred to as R-tipping. Another more recent study (11) has discovered how the stress response in bacteria is determined by the rate of environmental change. An increase in environmental stress leads to a single uniform pulse of alternative sigma factor σB activation, a general stress response pathway, with amplitude depending on the rate at which the stress increased. It is found that faster stresses lead to larger and sharper activation of σB, reflecting the fact that the activation process is rate-dependent. A question naturally arises how rate dependent signaling will affect the immune cell-fate selection via a differentiation of progenitor cells. We have studied these phenomena in the most paradigmatic switch responsible for the differentiation of immune cells, the GATA1-PU.1 switch. Moreover, we have considered how the shape of external signals may have an impact in decision-making process. The paper is structured as follows, we review the model of Huang et al. (2) and investigate, in addition to the symmetric scenario, the asymmetric scenario in two ways: (i) under the effect of asymmetric change of parameters; and (ii) under the effect of external signals, using two kinds of signals (see Materials and Methods). Furthermore, we will test the effect of parameter sweeping speed on the distribution of trajectories in the attractors of the dynamical system.
2. Materials and Methods
2.1. The GATA1-PU.1 Gene Regulatory Circuit
where X1 and X2 are the concentrations of two transcription factors GATA1 and PU.1, respectively. These equations model the dynamics of self-activation and cross-inhibition with Hill functions (12). The parameters a1, a2 represent self-activation rates, the parameters b1, b2 are basal expression rates, k1, k2 are deactivation rates, the parameters r’s are thresholds at which the inflection point in the Hill function occurs, and n is the Hill coefficient. The first terms of equations (1) and (2) give the contribution from self-activation, while the second terms measure the effect of cross-inhibition on basal activation rates, and the third terms the degradation. To take account of intrinsic gene expression stochasticity, we consider the differential equations (1) and (2) in the Langevin form by adding multiplicative noise terms (the last ones) where and stand for a Gaussian noise and depend on X1,2 as suggested in Ref. (13). These noise terms model the contribution of intrinsic noise which is unavoidable in biological systems. External cell signaling can be included in the model as follows
where S1 and S2 represent external signals to the genetic switch. Here, we are interested in two generic forms of signals:
• Linear signals: In this form (7) the external signals may have different rising times but they are equal in the steady state at Smax = 10 (see Figure 2A). For the sake of simplicity we assume that S1 reaches to the steady state faster than S2, and thus the rising time T1 of S1 is smaller than the rising time T2 of S2. They both increase linearly with time according to
The difference between S1 and S2 and the maximal difference A (Figure 2B) are defined as follows
• Adaptation form of signals: As suggested in Ref. (14) to achieve biochemical adaptation the signals have transient growth stage where they reach to their maxima, and decay stage where they decay and saturate to their steady states (see Figure 3). As for the first form, S1 has a rising time, θ1, smaller than S2, θ2, and the value of saturation is 10. They have the following form
where h1, h2 control the amplitude of signals, and θ1, θ2 are scale parameters. The second terms control the saturation of the signals to the value v = 10 (selected value). The maxima occur at tmax = θ1,2. Consequently, we have chosen θ (θ1 or θ2) to be the value which determines the speed since as we increase θ, which represents the rising time, we increase the time at which the maximum occurs. In other words, we decrease the speed of the signal variation.
Figure 1. GATA1-PU.1 genetic switch with and without external signals. (A) The isolated switch consists of two transcription factors GATA1 and PU.1 that activate themselves while inhibit each other’s expression. (B) The exposure of the same switch in (A) to two external signals S1 and S2.
Figure 2. Linear form of external signals in GATA1-PU.1 genetic switch. (A) Two external signals S1 and S2 with different rising times but equal steady states at Smax = 10. (B) The difference between the external signals with maximal asymmetry at A.
Figure 3. Adaptation form of external signals in GATA1-PU.1 genetic switch. The external signals S1 and S2 have different rising times but equal steady states at v = 10. Note that at the end of the signaling the system is structurally symmetric.
2.2. Testing of the Parameter Sweeping Speed
To test the effect of speed, we compute the ratio R numerically using
It represents the ratio between trajectories or cells which go to the top (or to the upper branch) of the bifurcation diagram, and trajectories that go to both upper and lower branches during simulation. Obviously, R = 1 if all cells choose the upper branch in the decision of their fate, R = 0.5 if the proportions of cells between two branches are equal, and R = 0 if all cells prefer the lower branch.
Heun’s method is used for solving the differential equations. In simulation of stochastic differential equations we have used Matlab, and all bifurcation diagrams and nullclines were generated in XPPAUT. δ(t) is an integration step size.
3.1. GATA1-PU.1 Genetic Switch without External Signals
This switch (Figure 1A) represents a paradigm for gene regulatory networks that govern the differentiation (2). It consists of two transcription factors GATA1 and PU.1 with self-stimulation and cross-inhibition. GATA1 is a master regulator of the erythroid lineage, and PU.1 is a master regulator of the myeloid lineage, and the two lineages arise from a common myeloid progenitor cell (1, 15).
3.1.1. Bifurcation analysis for symmetric scenario
In the symmetric scenario, the parameters of the model are changed symmetrically with respect to X1 and X2. Hence, the rates of self-activation, cross-inhibition, deactivation, and thresholds are equal for both transcription factors (see Materials and Methods). Then, this scenario is divided into two parts depending on the kind of bifurcation which results in during a change of the parameters.
• Supercritical pitchfork bifurcation: This type of bifurcation can occur when b is increased from 0.5 to 1 (Figure 4A), or when r is decreased from 1.8 to 1.2 (Figure 5C). In this kind of bifurcation, a transition occurs from monostability to bistability. The monostable state represents progenitor cell in undifferentiated state and has the ability to differentiate into two different fates. At this state, both transcription factors in the network are produced at approximately equal levels as it can be seen from the intersection point of nullclines in (Figure 4B). At the differentiation process, the progenitor cell is destabilized and two new attractors appear with equal basins of attraction (Figure 4C).
• Subcritical pitchfork bifurcation: This type of bifurcation occurs for many parameter changes. It can happen when k is changed from 1 to 1.5 (Figure 5A), when a is decreased from 1 to 0.5 (Figure 5B), when b is increased from 0.3 to 0.4 (Figure 5D), and when r is increased from 0.5 to 1 (Figure 5C). In this kind of bifurcation, a transition occurs from tristability to bistability (Figures 5E,F). In this situation, the progenitor cell (metastable state) coexists with the two fates, and the two transcription factors are expressed at equal or low levels. At the bifurcation process, it becomes unstable and makes discontinuous transition to either erythroid or myeloid fates with equal basins of attraction.
Figure 4. Supercritical pitchfork bifurcation diagrams for symmetric scenario. Bifurcation diagram (A) and nullclines at the beginning (B) and end (C) of the bifurcation. For all diagrams, n = 4, r = 0.5, a1 = a2 = 0.01, k1 = k2 = 1. The solid lines indicate stability, while dashed lines indicate unstable branches.
Figure 5. Subcritical pitchfork bifurcation diagrams for symmetric scenario. Bifurcation diagrams (A–D) and nullclines at the beginning (E) and end (F) of the bifurcation diagram (D). For all n = 4. For (A) a = 1, b = 1, r = 0.5, (B) b = 1, k = 1, r = 0.5, (C) a = 1, b = 1, k = 1, (D–F) a1 = a2 = 1, k1 = k2 = 1.5. In (C) there is also supercritical pitchfork bifurcation.
3.1.2. Bifurcation analysis for asymmetric scenario
Here, the parameters of the model are changed asymmetrically with respect to X1 and X2. For example, we can increase one of the parameters and keep the other constant, or decrease one of the parameters and keep the other constant, or both. This asymmetric change will cause symmetry-breaking in the bifurcation diagrams and makes one of the attractors more favorable than the other. Similarly, we note two types of bifurcation:
• Supercritical pitchfork bifurcation: In order to get this kind of bifurcation with symmetry-breaking, we increase a1 and k2 (see Materials and Methods for their definitions) and as a result, X1 is increased. Then, the bifurcation occurs when b is changed from 0.6 to 1 (Figures 6A–C). Now, the uncommitted progenitor cell represented by monostability is not in the middle but at the point where X1 is higher. After bifurcation, the erythroid fate becomes dominant since it has a larger basin of attraction to the right of the separatrix (Figure 6C).
• Subcritical pitchfork bifurcation: The bifurcation occurs when b is varied from 0.3 to 0.4. This gives imperfect subcritical pitchfork bifurcation (Figure 7A). The change in system behavior from tristability to bistability is depicted in (Figures 7B,C). At the progenitor cell, both transcription factors have low levels but the progenitor cell is not exactly in the middle. After bifurcation, one of the attractors corresponding to erythroid lineage becomes dominant as a result of increasing self-activation of GATA1.
Figure 6. Asymmetric scenario. Supercritical pitchfork bifurcation diagrams. Bifurcation diagram (A) and nullclines at the beginning (B) and end (C) of the bifurcation. The parameters are n = 4, r = 0.5, a1 = 0.2, a2 = 0.01, k1 = 1, k2 = 1.1.
Figure 7. Asymmetric scenario. Subcritical pitchfork bifurcation diagrams. Bifurcation diagram (A) and nullclines at the beginning (B) and end (C) of the bifurcation. The parameters are n = 4, r = 0.5, a1 = 1.2, a2 = 1, k1 = 1.5, k2 = 1.6.
3.1.3. Trajectories and the effect of parameter sweeping speed
To investigate the effect of the different speeds of the parameter sweeping we concentrate on the asymmetric supercritical pitchfork bifurcation, and similar results can be seen in the other kind of bifurcation. The graphical solutions of X1 and X2 after solving the differential equations (see equations (1) and (2) in Materials and Methods) are shown in (Figure 8A). As the time increases, the values of X1 increase and the values of X2 decrease. In fact, for small values of noise, this is the expected behavior from the dominance of the erythroid attractor.
Figure 8. Trajectories and parameter sweeping speed. (A) Time evolution of X1 and X2 in the asymmetric supercritical pitchfork bifurcation. (B) The effect of increasing the speed of crossing the critical region on the distribution of trajectories in the attractors for 10000 iterations. As the speed is increased, the ratio R changes from 1 to 0. Hence, increasing the speed causes a large switch from the favorable attractor to the other one. Parameters are a1 = 0.2, a2 = 0.01, k1 = 1, k2 = 1.1, n = 4, r = 0.5. Also, in (A) σ2 = 0.01, (B) σ2 = 0.5.
To examine the effect of the speed with which the system crosses the critical region, we vary b linearly with time according to b(t) = αt, where α is the slope, and compute the ratio R (see Materials and Methods). The result is shown in (Figure 8B). For low speeds, the ratio R is high which means that most of the cells choose the erythroid lineage due to the produced asymmetry, and this lineage leads to and include red blood cells. On the other hand, as we increase the speed, this ratio tends to zero. Two conclusions follow from this behavior. First, large speeds reduce the effect of asymmetry gradually and convert the favorable attractor completely when the ratio tends to zero. Second, R = 0 means that the myeloid fate becomes more favorable by cells. The myeloid fate leads to the immune cells of the immune system (16).
3.2. GATA1-PU.1 Genetic Switch Under External Signaling
To elucidate the effect of external signals on the dynamics of the switch, we consider external signals acting upon the switch (Figure 1B), see also equations (3) and (4). The external signals enhance the activation of X1 and X2. Figure 9 highlights the bifurcation in the parameter space (S1, S2) for the chosen set of parameters. The borders separate between the regions of monostability and the region of bistability, and this indicates to the existence of supercritical pitchfork bifurcation under the two following scenarios.
Figure 9. Two-parameter bifurcation diagram. The bifurcation in the parameter space (S1, S2), where S1 and S2 are external signals in the genetic switch. The borders separate between the regions of monostability I and the region of bistability II. Parameters are a1 = a2 = 0.05, b1 = b2 = 0.45, r = 0.5,k1 = k2 = 1.
3.2.1. Bifurcation analysis for symmetric scenario
Under this scenario, both signals S1 and S2 are equal. The nullclines in (Figures 10A,B) exhibit the bifurcation from monostability to bistability. This symmetry will give us near-symmetric bifurcation diagram (Figure 10C) with progenitor cell located in the middle and have equal probabilities to choose between erythroid (upper attractor) and myeloid (lower attractor) fates.
Figure 10. Nullclines and bifurcation diagrams with symmetric and asymmetric external signals. For all figures, a1 = a2 = 0.05, b1 = b2 = 0.45, r = 0.5, k1 = k2 = 1, n = 4. (A) Nullclines for S1 = S2 = 1 show one stable steady state, (B) nullclines for S1 = S2 = 4 show bistability, (C) near-symmetric supercritical pitchfork bifurcation, (D) nullclines for S1 = 3, S2 = 1 show one stable steady state shifted to the right, (E) nullclines for S1 = 6, S2 = 4 show bistability with larger basin of attraction to the right of the separatrix (almost diagonal line), (F) imperfect supercritical pitchfork bifurcation due to the asymmetry between the external signals.
3.2.2. Bifurcation analysis for asymmetric scenario
In contrary to the above scenario, now the signals have different parameters. As a result, the monostable state (Figure 10D) is at the point where X1 is higher since S1 which acts on X1 is larger. After bifurcation, the attractor at which X1 is high has a larger basin of attraction (Figure 10E). We can note in (Figure 10F) how this asymmetry produces symmetry-breaking in the bifurcation diagram and so the decision of the cell will be biased.
3.2.3. Trajectories and speed-dependent cellular decision making
To study how signal asymmetry, noise, and decision making will result in the dependence of the parameter sweeping speed we have considered with two kinds of signals (See Materials and Methods):
• Linear signals: The signals are shown in (Figure 2A). The asymmetry between the two signals is transient and the symmetry is retained after some time (Figure 2B). The behavior of trajectories of X1 and X2 under the influence of this form of signals is shown in (Figure 11A). As the time increases, the values of X1 increase and the values of X2 decrease. Hence, trajectories of X1 and X2 choose the attractor at which X1 is higher since S1 is faster. Next, to test the effect of increasing the speed on choosing the attractors (stable steady states) of genetic switch in the presence of noise and transient asymmetry A, we vary T1 in with constant values of A and Smax, and T2 will be changed according to the formula (7). With increasing the speed (Figure 11B), the ratio R tends to 0.5. Thus, increasing the speed increases the symmetry between erythroid and myeloid lineages and reduce the effect of asymmetry which is produced by external signals.
• Adaptation form of signals: The signals take the non-linear form shown in (Figure 3) and as for the linear form, S1 is faster than S2. The trajectories in this form behave almost like the first form (Figure 11C). To study the effect of the speed, we vary θ1 and θ2 such that θ1 is smaller than θ2. Then, we compute the ratio R and the result is depicted in (Figure 11D). It shows ratio tending to 0.5 as θ is increased. But increasing θ decreases the speed, so, surprisingly, now we regain the symmetry in the switch by decreasing the speed of external signals.
Figure 11. Trajectories and effect of the external signaling speed. In (A,C) the time evolution of X1 and X2 under the effect of linear and adaptation form of signals is shown, respectively. The values of X1 increase and the values of X2 decrease because S1 is chosen to be faster than S2. Hence, trajectories choose the attractor which has a larger value of X1. (B) The effect of increasing the speed of linear form of signals for 1000 iterations. As the speed is increased, the ratio R tends to 0.5. Thus, increasing the speed increases the symmetry in the switch. (D) The effect of speed with the adaptation form of signals. Decreasing the speed gives ratio R tending to 0.5. Surprisingly, now decreasing the speed increases the symmetry in the switch. Parameters in (A), (B) are A = 2.5, Smax = 10, (C,D) h1 = h2 = 10, v = 10, and for all we have a1 = a2 = 0.05, b1 = b2 = 0.45, r = 0.5, k1 = k2 = 1.
We have shown the importance of parameter sweeping speed when the gene regulatory circuit of immune cell differentiation is exposed to external factors that cause symmetry-breaking and make one of the attractors or fates more favorable than the other. In our study, symmetry-breaking is caused by three factors. The first factor is the asymmetric change of parameters which gives ratio tends to zero as the speed is increased (Figure 8B). This means we get large conversion from the favorite attractor, the erythroid lineage, to the myeloid lineage. The importance of this effect may appear in cases where the person has a problem with immunity due to the decrease in the production of immune cells, so even when there is a bias in the cell and this bias has the effect of choosing the attractor where GATA1 is upregulated, the cell can be forced to choose the attractor where PU.1 is upregulated by increasing the speed of crossing the critical region and so enhancing the production of immune cells. The second factor is linear form of signals (Figure 2A) and in this case we get ratio tends to 0.5 with increasing the speed (Figure 11B). This result may be important in situations that need symmetry between erythroid and myeloid cells, or when decreasing the probability of choosing the erythroid lineage is required. The third factor is represented by non-linear form of signals, i.e., signals describing biochemical adaptation (Figure 3). Here, decreasing the speed blinds the asymmetry and produce symmetry between the two lineages (Figure 11D). Taken together, the external signals, its shape, and its speed may have critical effects on choosing the attractors and affect the cell-fate determination.
Notably, we followed the model of Huang et al. (2) to study the differentiation into erythroid and myeloid fates. On the other hand, there is a scheme in Ref. (1, 17) that gives additional kinds of cells or lineages under each transcription factor. In this scheme, GATA1 is responsible for differentiation into erythroid or megakaryocyte cells, and PU.1 leads to either lymphoid lineage (B and T cells) which gives the Adaptive Immune Cells, or myeloid lineage (macrophages and granulocytes) that produces the Innate Immune Cells. So for this scheme, the importance of parameter sweeping speed is increased as the fate corresponding to high concentration of PU.1 is able to produce the different types of immune cells.
Of particular interest and agreement with our conclusions about the importance of external signals, Heuser et al. (6) have showed the crucial role of external signals in MN1 leukemia. They have investigated the requirement of FLT3 and c-Kit signals for MN1 leukemia. Overexpression of MN1 induces myeloid leukemia and blocks erythroid differentiation. FLT3 and c-Kit signaling direct MN1-expressing cells toward the myeloid lineage, so disruption of these signals may prevent leukemia. Interestingly, the disruption of these external signals doesn’t delay the disease latency but induces a switch from myeloid to erythroid lineage. Thus, the external signals can alter leukemia stem cell differentiation fates.
Many models have focused on the role of external signals in the differentiation process (3, 5) but they haven’t given any attention to the shape of signals or to the speed of these signals. Additionally, many works have made their studies limited to the symmetric scenario for the sake of simplicity (2, 18). But in this paper, we have studied the asymmetric scenarios and investigated the effect of external signaling speed on the system’s dynamics. As a prospect, it would be specially interesting to study the effect of speed on more complicated models and including other factors that may have a role in the differentiation process of hematopoietic stem cells, which can lead to better understanding of the immune system. Furthermore, an experimental evidence is needed to support the predictions from the mathematical models.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
This project was supported by the Deanship of Scientific Research (DSR), King Abdulaziz University (KAU), Jeddah, under grant No. (20/34/Gr). Alexey Zaikin acknowledges support from CR-UK and Eve Appeal (PROMISE-2016).
1. Laslo P, Pongubala JMR, Lancki DW, Singh H. Gene regulatory networks directing myeloid and lymphoid cell fates within the immune system. Semin Immunol (2008) 20:228–35. doi: 10.1016/j.smim.2008.08.003
3. Chickarmane V, Enver T, Peterson C. Computational modeling of the hematopoietic erythroid-myeloid switch reveals insight into cooperativity, priming, and irreversibility. PLoS Comput Biol (2009) 5(1):e1000268. doi:10.1371/journal.pcbi.1000268
4. Palani S, Sarkar CA. Positive receptor feedback during lineage commitment can generate ultrasensitivity to ligand and confer robustness to a bistable switch. Biophys J (2008) 95:1575–89. doi:10.1529/biophysj.107.120600
5. 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):e1000518. doi:10.1371/journal.pcbi.1000518
6. Heuser M, Park G, Moon Y, Berg T, Xiang P, Kuchenbauer F. Extrinsic signals determine myeloid-erythroid lineage switch in MN1 leukemia. Exp Hematol (2010) 38:174–9. doi:10.1016/j.exphem.2010.01.003
10. Ashwin P, Wieczorek S, Vitolo R, Cox P. Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system. Philos Trans A Math Phys Eng Sci (2012) 370:1166–84. doi:10.1098/rsta.2011.0306
17. Adolfsson J, Mansson R, Buza-Vidas N, Hultquist A, Liuba K, Jensen CT, et al. Identification of Flt3+ lympho-myeloid stem cells lacking erythro-megakaryocytic potential: a revised road map for adult blood lineage commitment. Cell (2005) 121:295–306. doi:10.1016/j.cell.2005.02.013
Keywords: GATA1-PU.1 switch, differentiation, immune cells, pluripotent cells
Citation: Alagha A and Zaikin A (2013) Asymmetry in erythroid-myeloid differentiation switch and the role of timing in a binary cell-fate decision. Front. Immunol. 4:426. doi: 10.3389/fimmu.2013.00426
Received: 28 August 2013; Paper pending published: 12 October 2013;
Accepted: 20 November 2013; Published online: 05 December 2013.
Edited by:Carmen Molina-Paris, University of Leeds, UK
Reviewed by:Koji Yasutomo, University of Tokushima, Japan
Maria L. Toribio, Spanish Research Council (CSIC), Spain
Copyright: © 2013 Alagha and Zaikin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Alexey Zaikin, Department of Mathematics and Institute for Women’s Health, University College London, Gower Street, London WC1E 6BT, UK e-mail: email@example.com