A Mathematical Model for the Reciprocal Differentiation of T Helper 17 Cells and Induced Regulatory T Cells

The reciprocal differentiation of T helper 17 (TH17) cells and induced regulatory T (iTreg) cells plays a critical role in both the pathogenesis and resolution of diverse human inflammatory diseases. Although initial studies suggested a stable commitment to either the TH17 or the iTreg lineage, recent results reveal remarkable plasticity and heterogeneity, reflected in the capacity of differentiated effectors cells to be reprogrammed among TH17 and iTreg lineages and the intriguing phenomenon that a group of naïve precursor CD4+ T cells can be programmed into phenotypically diverse populations by the same differentiation signal, transforming growth factor beta. To reconcile these observations, we have built a mathematical model of TH17/iTreg differentiation that exhibits four different stable steady states, governed by pitchfork bifurcations with certain degrees of broken symmetry. According to the model, a group of precursor cells with some small cell-to-cell variability can differentiate into phenotypically distinct subsets of cells, which exhibit distinct levels of the master transcription-factor regulators for the two T cell lineages. A dynamical control system with these properties is flexible enough to be steered down alternative pathways by polarizing signals, such as interleukin-6 and retinoic acid and it may be used by the immune system to generate functionally distinct effector cells in desired fractions in response to a range of differentiation signals. Additionally, the model suggests a quantitative explanation for the phenotype with high expression levels of both master regulators. This phenotype corresponds to a re-stabilized co-expressing state, appearing at a late stage of differentiation, rather than a bipotent precursor state observed under some other circumstances. Our simulations reconcile most published experimental observations and predict novel differentiation states as well as transitions among different phenotypes that have not yet been observed experimentally.


Introduction
CD4 + T cells are important components of the adaptive immune system in higher vertebrates. By producing various cytokines, they perform critical functions such as helping B cells to produce antibodies, activating CD8 + cytotoxic T cells, enhancing the innate immune system, and suppressing the immune response to avoid autoimmunity [1,2,3]. In peripheral tissues, such as lymph nodes, blood and sites of infection, antigen-inexperienced (naïve) CD4 + T cells can differentiate into effector cells of specialized phenotypes upon stimulation by cognate antigen delivered to the T cell receptor by Antigen Presenting Cells (APCs). Proliferation and differentiation of activated naïve T cells depends on their particular cytokine microenvironment. These specialized effector T cells produce distinct cytokine profiles tailored for their specialized functions. Also, they express lineage-defining transcription factors (''master regulators''). In general, high expression level of a particular master regulator is observed only in cells of a particular lineage, and the overexpression of a particular master regulator induces the production of the corresponding lineage-defining cytokines [4,5].
The fate of a naïve CD4 + T cell was traditionally thought to be either T helper 1 (T H 1) cell or T helper 2 (T H 2) cell [6]. In the last decade, a third type of T helper cell (T H 17), derived from naïve CD4 + T cells, was discovered [7]. T H 17 cells produce interleukin-17A (IL-17A), IL-17F and IL-22 as their lineage-defining cytokines, and the retinoic acid receptor-related orphan receptor gamma t (RORct) transcription factor is considered the master regulator of this lineage [8,9]. In addition, naïve CD4 + T cells were found to be able to differentiate into a fourth lineage of (regulatory) T cells, which were called induced regulatory T (iT reg ) cells to distinguish them from natural regulatory T (nT reg ) cells, which differentiate in the thymus instead of the periphery [10]. iTreg cells are characterized by producing IL-10 and transforming growth factor-b (TGF-b) and highly expressing forkhead box P3 (Foxp3) transcription factor as their master regulator [11]. T H 17 cells are pro-inflammatory because they secret cytokines that promote inflammation, whereas iT reg cells are anti-inflammatory because their lineage-defining cytokines can reduce the inflammatory response.
The differentiation pathways of naïve T cells into T H 17 and iT reg lineages are closely related. First, stimulation by TGF-b is necessary for the differentiation of both lineages [12]. The differentiation of T H 17 and iT reg cells are reciprocally regulated in the presence of TGF-b, i.e. inhibiting the differentiation pathway of one lineage will result in activation of the pathway for the other lineage. This is due to the mutual antagonism between RORct and Foxp3. Furthermore, polarizing signals, such as IL-6 and retinoic acid, can induce the differentiation of one lineage and repress that of the other one [12]. Nonetheless, differentiated iT reg cells can be reprogrammed into T H 17 cells in an appropriate cytokine environment [13], suggesting significant plasticity of these two lineages. In addition, stable co-expression of their master regulators (RORct and Foxp3) is observed both in vivo and in vitro [14,15]. Interestingly, these double-expressing cells were found to possess either regulatory or dual (regulatory and proinflammatory) functions in vivo [14,15].
Perhaps the most intriguing phenomenon is that antigenactivated naïve CD4 + T cells treated with TGF-b alone give rise to a heterogeneous population, which may include three phenotypes (Foxp3-only, RORct-only, and double-expressing cells) at an intermediate TGF-b concentration [16], or two phenotypes (RORct-only and double-expressing cells) at a higher TGF-b concentration [15]. In combination with TGF-b, IL-6 can induce the differentiation of RORct expressing cells, whereas all-trans retinoic acid (ATRA) can induce the differentiation of Foxp3 expressing cells [16,17] (Figure 1). All of these in vitro derived phenotypes can be observed in vivo, and at least some of their respective functions have been demonstrated, suggesting that these in vitro differentiation assays provide important clues to our understanding of the development of T H 17 and iT reg cells in the body.
Mathematical modeling has contributed to our understanding of the differentiation of T H 1 and T H 2 cells [18,19,20,21,22,23,24]. Höfer et al. first demonstrated that the dynamics of the key transcription factors can govern the robustness of the lineage choice and maintenance [18,19]. Yates et al. later related transcription factor dynamics to the mix of T H 1 and T H 2 cells in a population of differentiating T cells [20]. Recently, Bonneau et al. [25] have proposed a Boolean-network model of the comprehensive repertoire of CD4 + T cell phenotypes, including T H 17 and iT reg cells. Drawing inspiration from these earlier models, we have sought to explain, with a computational model, the remarkable heterogeneity of the T H 17-iT reg reciprocaldifferentiation system.
In terms of this model, we show that a population of naïve CD4 + T cells, with some small cell-to-cell variability, can differentiate into a heterogeneous population of effector cells with distinct phenotypes upon treatment with the primary differentiation signal (TGF-b). Polarizing signals, such as IL-6 and ATRA, can skew the differentiation to one or two phenotypes. A control system with these properties can generate functional diversity of the induced cell populations and can be regulated with great flexibility by diverse environmental cues. In addition, the model suggests how treatment with different concentrations of TGF-b may favor different responding phenotypes, and how conversions among these phenotypes may be guided. Finally, the model gives a new quantitative explanation for double-expressing cells, suggesting that they are 're-stabilized co-expressing' cells rather than transient intermediate cells in the differentiation pathway. The model predicts that double-expressing cells should appear at a relatively late stage of the differentiation process, and they may be intended for specific functions. In all, our model provides a novel mathematical framework for understanding this reciprocal differentiation system, and it gives new insights into the regulatory mechanisms that underlie the molecular control of certain immune responses.

Results
A model with symmetrical interactions predicts three differentiated phenotypes of CD4 + T cells induced by TGF-b To illustrate our basic idea, we first construct a model of a simple and perfectly symmetrical regulatory network (Figure 2A). In the Methods section we describe how this network is converted into a pair of nonlinear ordinary differential equations (ODEs) for

Author Summary
In order to perform complex functions upon pathogenic challenges, the immune system needs to efficiently deploy a repertoire of specialized cells by inducing the differentiation of precursor cells into effector cells. In a critical process of the adaptive immune system, one common type of precursor cell can give rise to both T helper 17 cells and regulatory T cells, which have distinct phenotypes and functions. Recent discoveries have revealed a certain heterogeneity in this reciprocal differentiation system. In particular, treating precursor cells with a single differentiation signal can result in a remarkably diverse population. An understanding of such variable responses is limited by a lack of quantitative models. Our mathematical model of this cell differentiation system reveals how the control system generates phenotypic diversity and how its final state can be regulated by various signals. The model suggests a new quantitative explanation for the scenario in which the master regulators of two different T cell lineages can be highly expressed in a single cell. The model provides a new framework for understanding the dynamic properties of this type of regulatory network and the mechanisms that help to maintain a balance of effector cells during the inflammatory response to infection.
the time rates of change of Foxp3 and RORct. The rate functions for this model contain 12 kinetic parameters, whose basal values are specified in the Methods section (Table 1) for the ''symmetrical model without intermediates''. The solution of these ODEs for the basal values, and with [TGF-b] = 0, evolves to a stable steady state where both RORct and Foxp3 have a low level of expression (RORct low Foxp3 low ). This steady state corresponds to a naïve CD4 + T cell ( Figure 3A). In the presence of a sufficient TGF-b signal, the regulatory network might evolve to one of three other steady states, namely RORct high Foxp3 low , RORct low Foxp3 high and RORct high Foxp3 high states, corresponding to RORct-only, Foxp3-only and double-expressing phenotypes. Note that these stable steady states are also referred as 'cell fate attractors' in some other studies, and this concept facilitates our understanding of cell lineage choice and reprogramming (reviewed in [26]). Figure 3B shows a scenario in which the TGF-b signal triggers the formation of a tri-stable system. In this particular case, the RORct low Fox-p3 low state is no longer a stable steady state, and naïve cell, which was previously stabilized in the RORct low Foxp3 low state, would differentiate into the RORct high Foxp3 high state, whose basin of attraction (the white region in Figure 2B) contains the naïve state of the cell.
However, cell-to-cell variability can produce other results. We interpret cell-to-cell-variability as small deviations of the parameter values from their basal settings in Table 1. The basal settings correspond to the behavior of an ''average'' cell, but any particular cell will deviate somewhat from this average behavior. As consequences of the changing parameter values in any particular cell, the position of the RORct low Foxp3 low state changes, the boundaries of the basins of attractions change, and the fate of the naïve cell may change. The naïve T cell will differentiate into the stable steady state in whose basin of attraction it lies. That is, depending on the precise parameter values of the cell, its RORct low Foxp3 low state may lie in any of the three basins of attraction of the TGF-b-stimulated system. Figure 3C depicts three cells in the population that adopt three different fates because of the variability among them. With a random sample of cells, each of the three differentiated states can be populated by a significant fraction of cells ( Figure 3D). Although cell-to-cell variability does not make large changes in the position of the RORct low Foxp3 low state, it has a dramatic influence on the basins of attraction of the stable steady states, which determines the fate of the cell once the differentiation signal is turned on.
Since the system has four distinct steady states that correspond to four distinct phenotypes, we next looked for the relationships among these steady states using bifurcation analysis of an average cell. Because of the symmetrical nature of the interactions, an average cell exhibits sub-critical pitchfork bifurcations with TGF-b concentration as the control parameter ( Figure 4A). (The notion of a pitchfork bifurcation was used earlier, in references [27,28], to explain a system of hematopoietic cell differentiation in which multiple lineages might be adopted.) Notably, the RORct low Fox-p3 low state is only stable at low TGF-b concentration. At an intermediate concentration of TGF-b (,0.25 units in Figure 3A), the system bifurcates into two lineage-specific branches, corresponding to RORct high Foxp3 low and RORct low Foxp3 high states. The fourth type of stable steady state (RORct high Foxp3 high ) appears at higher TGF-b signal strength (.0.37 in Figure 3A  transcription factors Foxp3 and RORct ( Figure 2B). In this case, the system can be tri-stable even at high concentrations of TGF-b, and the total conversion of single-expressing cells into double-expressing cells would not occur. Instead, co-existence of the three phenotypes in comparable fractions might be observed over a wide range of [TGF-b] ( Figure 4B).

A model with asymmetrical interactions provides a better account of the regulatory functions of TGF-b during the coupled differentiation of T H 17 and iT reg cells
We next considered an asymmetrical model in which the network topology and parameter values differ from the symmetrical model. In the model with perfect symmetry, we assumed that the inhibitions between Foxp3 and RORct are equally strong, which is not supported by existing experimental evidence. In fact, Foxp3 is better known for its inhibitory function on IL-17, a downstream effector of RORct, as demonstrated by Williams and Rudensky [29]. Therefore, we revised our model by removing the direct inhibition of RORct expression by Foxp3 and adding the inhibition of IL-17 expression by Foxp3. This revised model, with broken symmetry ( Figure 1C, Table 1  phenotypes in comparable fractions have been observed [16]. In addition, the maximum percentage of Foxp3 single-positive cells was observed at some lower concentration of TGF-b. As [TGF-b] was increased, the percentage of Foxp3 single-positive cells decreased, accompanied by a concordant rise in the percentage of RORct-expressing cells [16]. At higher concentrations of TGFb, RORct-only cells and double-expressing cells were found to coexist in comparable percentages [15].
Our model not only validates existing published data on the coexistence of two or more phenotypes in mixed T helper cell populations but also predicts that increasing TGF-b concentration will cause the transformation of Foxp3 single-positive cells into RORct-expressing cells. Conversely, decreasing TGF-b concentration might result in the reverse transformation.

Our model accommodates the observed effect of IL-6 skewing T cells into a 'RORct-only' phenotype
We next simulated the influence of IL-6 on this reciprocal differentiation system. In the asymmetrical model ( Figure 3C), IL-6 activates STAT3, which favors production of RORct over Foxp3. In this model, IL-6 will not trigger differentiation in the absence of TGF-b. However, IL-6 significantly increases the fraction of RORct-only cells over a wide range of TGF-b concentrations ( Figure 4A). Also, it stimulates some of the cells in the (simulated) population to produce IL-17. These results are consistent with the observations of a few groups [13,16]. In particular, Zhou et al. observed that low level TGF-b favors the RORct-only phenotype and IL-17 production, whereas higher concentrations of TGF-b inhibit the production of IL-17. They also reported that the decrease of IL-17 production at higher TGF-b concentration is accompanied by an increase of Foxp3expressing cells. We see this phenomenon in our simulation, and we further suggest that the decrease of RORct-only cells, or the increase of the double-expressing cells, accounts for the reduced production of IL-17 at high TGF-b concentration, because double-expressing cells are known to be much less effective in producing IL-17 than the RORct-only cells, at least in this type of in vitro assay with TGF-b and IL-6 [15,16]. However, Zhou et al. observed a pronounced inhibition of IL-17 production at higher TGF-b concentration even when Foxp3 expression had not been remarkably raised [16]. This discrepancy suggests that high TGFb level may trigger Foxp3-independent repression of IL-17 production.
Both the observations by Zhou et al. and our simulations demonstrate that only a minor fraction of RORct-only cells exhibit IL-17 production even in the presence of IL-6. In fact, this is not an idiosyncratic phenomenon. Mariani et al. recently discovered that only a subset of T H 2 cells produce IL-4 due to cellto-cell variability [30], suggesting that the production of lineagespecific cytokines in T helper cells can be controlled by stochastic mechanisms.

Our model accommodates the effect of ATRA skewing T cells into a Foxp3-expressing phenotype
In the asymmetrical model ( Figure 3C), ATRA favors production of Foxp3 over RORct. Hence, in our simulation of TGF-b+ATRA stimulation, we found that the percentage of Foxp3-only cells and double-expressing cells significantly increased as compared to TGFb alone (compare Figure 4B to Figure 3C). Like IL-6, ATRA did not trigger differentiation by itself. We next checked if ATRA can suppress the polarizing effect of IL-6. In our simulation, ATRA was effective in reducing the IL-6 induced production of IL-17. In addition, at high TGF-b concentration, ATRA significantly decreased the percentage of RORct-only cells, and resulted in a population with comparable fractions of RORct-only cells and double-expressing cells ( Figure 5C). All of these simulation results are consistent with published data [13,15,17,31]. Our model suggests that ATRA can significantly increase the percentage of Our model predicts that IL-6 may reprogram iT reg cells to IL-17 producing cells, while ATRA may prevent this reprogramming effect With our model, we next checked whether IL-6 could reprogram differentiated iT reg cells into T H 17 cells. We first induced a population of naïve CD4 + T cells to differentiate into a population dominated by 'Foxp3-only' cells with an intermediate level of TGF-b (0.28 units). After the cells came to their Foxp3only steady state, we raised the IL-6 signal to 10 units and continued the simulation. We found that almost all the cells expressing Foxp3 before adding IL-6 stopped producing Foxp3 upon the treatment with IL-6, and a subset of 'RORct-only' cells dominated the population. A fraction of these RORct-only cells produced IL-17 ( Figure 6A, left panel).
When we induced the differentiation of iT reg cells with TGF-b+ATRA and performed the same reprogramming simulation, we found that ATRA did not prevent the repression of Foxp3 expression by IL-6 significantly. However, ATRA prevented the formation of IL-17 producing cells ( Figure 6A, right panel). The reprogramming capability of IL-6 and the inhibitory effect of ATRA have been observed by Yang et al. [13].
Analyzing the concentration dependence of these reprogramming effects, we found that a high level of IL-6 may exclusively down-regulate Foxp3 expression ( Figure 6B, left panel) whereas a high level of ATRA may predominantly prevent IL-17 expression ( Figure 6B, right panel). Interestingly, when both of these factors are present in high concentration, our model predicts that, although most cells exhibit high expression of RORct, there are almost no IL-17-producing cells in the population. Future experimental studies are warranted to confirm these intriguing predictions. Table 2 summarizes the observations that are in agreement with our simulation results and the testable predictions that we have made based on the bifurcation analyses and signal-response curves.

Discussion
Previous mathematical models have shown how differentiation signals can trigger a robust switch during the development of T H 1 or T H 2 cells [18,19,20,21,22,23,24]. In particular, earlier modeling studies by Höfer et al. demonstrated how the interactions among transcription factors can create a memory for T H 2 lineage commitment and govern the choice of T H 1 and T H 2 lineages [18,19]. These studies focused on the dynamics of transcription factors within a single (average) cell, but the authors also pointed out that cell-to-cell variability in a CD4 + T cell population can be modeled mathematically by introducing parametric variations to the ordinary differential equations (ODEs). In addition to modeling molecular interactions, the study by Yates et al. related the dynamics of transcription factors to the phenotypic composition of T H 1 and T H 2 cell populations [20]. The authors built comprehensive ODE-based models which take into account cell proliferation, intercellular communication, and cell-to-cell variability. Yates et al. modeled cell-to-cell variability by variations in initial conditions, but we consider parametric variations to be a more important source of cell-to-cell variability (see Methods).
The reciprocal differentiation of T H 17 and iT reg cells, although a relatively new research field, has already been shown to exhibit many interesting and unique features, and yet it has not been studied in quantitative detail using mathematical models. The work presented here reveals some of the intriguing regulatory mechanisms of this differentiation system. We showed that the four phenotypes of cells, corresponding to four different steady states of the dynamical system, are derived from a pitchfork bifurcation with certain degree of broken symmetry. A single primary differentiation signal, TGF-b, can give rise to multiple cell types with distinct functions, while other polarizing differentiation signals, such as IL-6 as ATRA, skew the system to particular type(s) of cells. If we regard TGF-b as tossing dice for the naïve cells, those polarizing signals may load the dice, although they may not toss the dice themselves. The remarkable advantage of this system is that functionally synergic cells could be generated simultaneously in desired fractions with some simple differentiation inducers.
Our model suggests that the double-expressing phenotype is a re-stabilized co-expressing state, which should be observed in relatively late stages of cell differentiation. Previously, van den Ham and de Boer found this type of state in a similar dynamical system, although they chose parameter values to avoid this state for their system [24]. With perfectly symmetrical models, some other groups described a double-expressing state as an intermediate state before the decision making switch, corresponding to some bipotent precursor cells [27,32,33]. For the T H 17-iT reg paradigm, it is also possible that these double-expressing cells are at an intermediate state that should be converted into singleexpressing cells at a later stage of the differentiation process. However, we do not favor this view for the following reasons. 1) A few studies have shown that the double-expressing cells are effective in repressing effector cell growth and/or secreting proinflammatory and anti-inflammatory cytokines [15,34]. It is not likely that a differentiation intermediate would perform any conspicuous function in the immune system. 2) There are a few reports demonstrating the conversion from iT reg cells to double-expressing cells [13,14], or from RORct-only cells to doubleexpressing cells [15], and to our knowledge it is not yet established that observable double-expressing cells can be converted into single-expressing cells. Assuming that differentiation from early stage to late stage is more readily to be observed than the 'dedifferentiation' process, these results indicate that the doubleexpressing cells might be at a differentiation stage later than the single-expressing states. 3) As shown in this report, there is a mathematical basis to support the double-expressing state appearing only at relatively high TGF-b concentration and some late differentiation stage, and the model is in accord with most published experimental observations. In addition, we are aware that the double-expressing cells are also observed for iT reg -T H 1 and iT reg -T H 2 paradigms [3]. Therefore, the framework presented here may be helpful for understanding iT reg cells that express T- bet or GATA3 as well. Interestingly, conversion of Foxp3expressing iT reg cells to Foxp3/T-bet double-expressing cells has been reported [35]. In fact, these double-expressing cells may play very specific and indispensable roles in controlling inflammation. Chaudhry et al. have found that iT reg cells require STAT3 for their suppressive function on T H 17, and not on other lineages [36]. Koch et al. discovered that the T-bet expression is required for the function of iT reg cells during T H 1-mediated inflammation [35]. These results suggest that there are subpopulations of iT reg cells expressing various master regulators of T helper cells, and they are tailored for different functions [3]. Therefore, the doubleexpressing cells might be terminally differentiated effectors performing specific suppressive functions. It is possible that the Foxp3-only cells, which mainly appear at low TGF-b concentration, could serve as precursors or reservoir for different terminal effectors, in addition to their general suppressive functions.
Although the detailed physiological significance of this delicate differentiation system is yet to be discovered, Lochner et al. have already demonstrated in mice that, during infections and inflammation, the number of IL-17 producing RORct + cells and double-expressing cells increased in remarkably comparable proportions [15]. This suggests the need for balance between different cell types in response to pathogenic challenges. A single differentiation network that gives rise to multiple phenotypes might be crucial for the maintenance of such balance. Furthermore, it is worth highlighting the common features shared by the T H 17-iT reg differentiation system and the differentiation control systems of hematopoietic cells and of stem cells [27,28,37]. Functionally, these systems have the potential to generate multiple phenotypes in a single differentiation event, and these phenotypes may play synergic roles under certain physiological conditions. In addition, it has been shown that cell-to-cell variability within clonal populations makes significant contributions to the stochasticity of lineage choice in stem cells [38]. This is also concordant with our basic assumptions.
Pitchfork bifurcations (with broken symmetry) may be the underlying mechanism generating variable phenotypes in these dynamical control systems. We will not be surprised if other cell differentiation systems possess similar properties. Recently, Heinz et al discovered that the 'priming factor' PU.1, which is required for both macrophage and B cell differentiation, is responsible for creating some of the lineage specific epigenetic markers by itself [39]. Therefore, it is possible that these priming factors not only drive the differentiation event, but also help to create a heterogeneous population of cells.
One limitation of our model is the assumption that the high concentration of TGF-b used by Lochner et al. is above the saturation concentration for TGF-b signaling [15]. We are cautious about extrapolating our model to even higher TGF-b concentration because there is no available experimental result for us to compare with. In fact, it is possible that at even higher TGFb concentration either the RORct-only phenotype or the doubleexpressing phenotype dominates the population, and the conver-  [38]. Nonetheless, when more experimental results become available, we should be able to pinpoint the missing pieces in this reciprocal differentiation system and make the mathematical model more helpful for our understanding of the system in detail.
Another limitation of this study is that we have neglected the effects of intercellular communication on the differentiation of CD4 + T cells. Cytokines secreted by T H 1 and T H 2 cells are known to influence the differentiation of neighboring T cells [40], and previous modeling work has highlighted the importance of these paracrine signaling effects [20]. Relevant to our work, the cytokines secreted by T H 17 and iT reg cells can influence the differentiation of a population of T cells, and this influence might

Low-intermediate-high Higher concentration of TGF-b inhibits IL-17 production
Observed in more extent [16] Inducing differentiation from naïve CD4 + T cells with TGF-b and ATRA Intermediate More Foxp3 expressing cells than with TGF-b alone Observed [17] Intermediate Foxp3-only phenotype is the major differentiated phenotype Prediction High Double-expressing phenotype is the major differentiated phenotype Prediction Inducing differentiation from naïve CD4 + T cells with TGF-b, IL-6 and ATRA High RORct-only and double-expressing phenotypes in comparable fractions. IL-17 production is much lower than with TGF-b and IL-6 Observed [15] Inducing differentiation from naïve CD4 + T cells to iT reg cells with TGF-b, and reprogramming the differentiated iT reg cells with IL-6 Intermediate Foxp3 expressing cells are reduced, and IL-17 producing cells appear in significant fraction.
Observed [13] Inducing differentiation from naïve CD4 + T cells to iT reg cells with TGF-b and ATRA, and reprogramming the iT reg cells with IL-6 Intermediate Foxp3 expressing cells are reduced, and no significant number of IL-17 producing cells can be observed.
Observed [13] Intermediate Most cells are in 'poised' state at which RORct expression is high, but no IL-17 is produced. be reflected in changes of the proportions of induced phenotypes. For example, both T H 17 and iT reg cells can produce TGF-b [41,42], which may increase the percentage of both type of cells, or induce the transition from single-expressing cells to doubleexpressing cells, and this may be causative for the transition observed by Lochner et al. [15]. However, it is not yet clear how important are paracrine signals via secreted cytokines compared to exogenous cytokine signals, with respect to T H 17 and iT reg differentiation. We leave the consideration of these factors for future work. In summary, we presented a novel mathematical model of T H 17-iT reg differentiation. Based on the model, we show how TGF-b can trigger the differentiation of naïve CD4 + T cells into a heterogeneous population containing RORct-only, Foxp3-only and double-expressing cells, and how polarizing signals can skew the differentiation to particular phenotype(s). The model suggests how the conversions among different phenotypes can be guided. Additionally, the model gives a new quantitative explanation for the double-expressing cells, which should appear only at a late stage of the differentiation process. Our model provides new insights into the regulatory mechanisms that underlie the molecular control of certain immune responses.

Methods
We constructed our mathematical model based on known interactions among key molecules in the differentiation system of T H 17 and iT reg cells. For illustrative purposes, we first consider a 'symmetrical' model in which the lineages of T H 17 and iT reg have identical corresponding interaction types and strengths. Then we added two intermediate proteins for transmitting TGF-b signals in this symmetrical model. Next, we modified our model so that it became asymmetrical, and we incorporated two other input signals. Using this last model, we compared our simulation results with some published experimental data and made several testable predictions.
In the symmetrical model ( Figure 2A) TGF-b upregulates both RORct and Foxp3, which has been demonstrated in a few published experiments [13,43]. The model includes the 'autoactivation' of both master regulators. Although there is no evidence for direct autoactivation of RORct and Foxp3, these relationships in our model represent known positive feedback loops in their respective pathways. One origin of these positive feedback loops is the epigenetic modifications observed in the promoter regions of RORct and Foxp3 in their respective lineages [44,45]. These epigenetic modifications recruit additional chromatin remodeling complexes that further stabilize those modifications and help to maintain the gene expression, thus forming positive feedback loops [46]. Additionally, master regulators can enhance their own production by autocrine effects. For example, RORct can induce production of IL-21 and IL-23 which further stimulate the expression of RORct, as suggested by Murphy and Stokinger [47]. The symmetric model also includes the cross-inhibition interactions between Foxp3 and RORct. Inhibition of Foxp3 by RORct is supported by the recent discovery that RORct acts as a transcriptional repressor of Foxp3 by binding to its promoter [48]. Although a few reports suggest a functional inhibition of RORct by Foxp3 [13,16,49], the presence of Foxp3 was shown to have no pronounced effect on the expression of RORct [50]. Our symmetrical model includes the inhibition of RORct by Foxp3, but we relaxed this assumption in our model with broken symmetry.
In the first version of our symmetrical model, TGF-b directly activates RORct and Foxp3. In the second version, we added intermediate proteins between TGF-b and the master regulators. It is known that Smad2, Smad3 and Smad4 mediate the TGF-binduced upregulation of Foxp3 [51,52], but the Smad proteins are dispensable for upregulation of RORct. It is still unclear how the TGF-b signal is transmitted to RORct [52]. Thus, in Figure 1B The model with broken symmetry also includes IL-17, which is activated by RORct and STAT3, and deactivated by Foxp3 and ATRA [8,13,16,29,53]. As a polarizing signal, IL-6 stimulates RORct and IL-17 production, and represses Foxp3 expression through the STAT3 pathway [54]. Conversely, ATRA upregulates Foxp3, downregulates RORct, and inhibits IL-17 production [17,31]. These relations are all included in our model with broken symmetry ( Figure 2C).
To model the T H 17-iT reg reciprocal-differentiation system, we use a generic form of ordinary differential equations (ODEs) that describe both gene expression and protein interaction networks [55,56,57]. Each ODE in our model has the form: i~1,:::,N X i is the activity or concentration of protein i. X i (t) changes on a time scale = 1/c i . X i (t) relaxes toward a value determined by the sigmoidal function, F, which has a steepness set by s i . The basal value of F, in the absence of any influencing factors, is determined by v o i . The coefficients v j?i determine the influence of protein j on protein i. N is the total number of proteins in the network. For example, the pair of ODEs for the first symmetrical model are: All variables and parameters are dimensionless. One time unit in our simulations corresponds to approximately 1 hour. All simulations and bifurcation analyses were performed with PyDSTool, a software environment for dynamical systems [58]. In the Supplementary Information we provide a Python module file (Text S1) for PyDSTool that completely defines the ODEs we are solving in each case, and an example script (Text S2) to reproduce bifurcation diagrams shown in Figure 4A.
All the experimental results to which our model has been compared were obtained with differentiation assays that lasted 2-5 days, and these results are essentially consistent from one experiment to another. Thus, we assumed that the observed, differentiated cell phenotypes after 2-5 days are representative of stable steady states in our model.
We have chosen to use generic (phenomenological) ODEs instead of a more detailed kinetic model of the biochemical reaction network because we lack sufficient mechanistic and kinetic information on the molecular interactions in the T H 17-iT reg reciprocal-differentiation system. To build a detailed biochemical model, based on mass-action or Michaelis-Menten kinetics, would require us to make many assumptions on the underlying mechanism and rate constants with little or no experimental evidence to back up these assumptions. In such a case, a phenomenological model seems more appropriate to us. A similar approach has been adopted in earlier theoretical studies of T cell differentiation by Mendoza and Xenarios [22], who used a sigmoidal function similar to our F(sW), and by van den Ham and de Boer [21], who used Hill functions in place of our F(sW). To be sure that our results are not overly dependent on our mathematical approach, we have re-formulated our 'symmetrical model without intermediates' using Hill functions and confirmed that the model exhibits four types of stable steady states as [TGFb] is varied. The basic features of the bifurcation diagrams and signal-response curves are similar, regardless of which formalism is used (details available upon request).
To account for cell-to-cell variability in a population, we made many simulations of the system of ODEs, each time with a slightly different choice of parameter values, to represent slight differences from cell to cell. We assumed that the value of each parameter conforms to a normal distribution with CV = 0.05 (CV = coefficient of variation = standard deviation/mean). The mean value that we specified for each parameter distribution is also referred as the 'basal' value of that parameter (see Table 1). In our bifurcation analysis of the dynamical system, we consider an imaginary cell that adopts the basal value for each of its parameters, and we define this cell as the 'average' cell. Note that none of the cells in our simulated population is likely to be this average cell, because every parameter value is likely to deviate a little (CV = 5%) from the basal value. Note, in addition, that our simulations sample a volume of parameter space around the 'average' cell, thereby probing the sensitivity/robustness of the differentiation process. Because we are varying all parameters simultaneously and randomly, this procedure is more indicative of robust behavior than standard sensitivity analysis, which involves estimating the partial derivative of some output property (e.g., steady state level of Foxp3) with respect to each parameter separately.
In order to simulate the induced differentiation process, we first solved the ODEs numerically with some small initial values of [RORct] and [Foxp3] state and with [TGF-b] = 0 (and, if applicable, other input signals, e.g. IL-6 and ATRA, = 0 as well). After a short period of time, each simulated cell will find its own, stable RORct low Foxp3 low steady state, corresponding to a naïve CD4 + T cell. Next, we changed [TGF-b] (and other input signals, if applicable) to a certain positive value and continued the numerical simulation. By the end of the simulation, each cell arrives at its corresponding 'induced' phenotype, which might vary from cell to cell because of the parametric variability of the population. To simulate the reprogramming effect, the concentration of IL-6 was raised after the cells were stabilized in the differentiated state. We made the simple definition that a protein is expressed when its level is greater than 0.5 units.
To check the effect of TGF-b concentration on the induced phenotypes, we ran a series of simulations for a group of 1000 cells with various values of [TGF-b] and plotted the percentages of cells that adopt each terminal phenotype, in order to generate a 'signalresponse' curve for a population of cells. Note that this signalresponse curve could only represent a series of induced differentiation experiments with various TGF-b concentrations instead of a single experiment with increasing concentration of TGF-b.
Our simulations of cell-to-cell variability are based on the assumptions that each cell follows a deterministic trajectory but that cells differ from one another in the precise values of the kinetic parameters that govern the deterministic trajectory. A similar approach was adopted by Höfer et al. in their model of transcriptional regulation of T lymphocytes [18]. An alternative view of stochasticity assumes that all cells are identical in terms of kinetic constants but they follow unique stochastic trajectories because of random fluctuations in the numbers of molecules of the dynamic variables. The truth is most likely a combination of these effects (parameter variation and molecular fluctuations), but we have adopted the parameter-variation approach for several reasons. First of all, we lack the sort of molecular details (e.g., the numbers of molecules of regulatory species per cell) needed for accurate stochastic simulations of molecular fluctuations. Second, it is unlikely that T cells are identical with respect to parameter values, and there is experimental evidence to the contrary. Peripheral naïve T cells undergo a complex developmental process in the thymus, where they likely inherit many stable cell-to-cell differences, possibly because of the great diversity of T cell receptor specificities generated by VJ or V(D)J recombination. Experiments on T cell differentiation are done by selecting cells with some common characteristics, but they may nonetheless differ in many other respects. Even monoclonal populations of mammalian cells (derived from a single progenitor cell) exhibit a distribution of properties that can affect cell fate determination [38]. Nonetheless, to be sure that our results are not overly dependent on our view of cell-to-cell variability, we have reformulated our 'symmetrical model without intermediates' as a pair of stochastic differential equations with additive white noise and confirmed that the SDEs generate signal-response curves similar to our results in Fig. 4A, bottom panel (details available upon request).
It is also reasonable to attribute variability among cells to different initial conditions for each simulation of the governing ODEs, as suggested by Yates et al. [20]. Since variations of initial conditions can also bias cells toward different phenotypes, we presume that this strategy will produce results similar to our own.

Supporting Information
Text S1 A module file that defines the ODEs for the three models. This is a Python module file that specifies the equations and the parameter values for the three models discussed in the paper. They can be used as inputs for simulations and analyses with PyDSTool. (TXT) Text S2 An example script for generating bifurcation diagrams. This is a Python script file that produces the 1parameter bifurcation diagram shown in Figure 4A. It requires PyDSTool and the module file that defines the ODEs (Text S1). (TXT)