Abstract
Background
CD4^{+} T cells have several subsets of functional phenotypes, which play critical yet diverse roles in the immune system. Pathogendriven differentiation of these subsets of cells is often heterogeneous in terms of the induced phenotypic diversity. In vitro recapitulation of heterogeneous differentiation under homogeneous experimental conditions indicates some highly regulated mechanisms by which multiple phenotypes of CD4^{+} T cells can be generated from a single population of naïve CD4^{+} T cells. Therefore, conceptual understanding of induced heterogeneous differentiation will shed light on the mechanisms controlling the response of populations of CD4^{+} T cells under physiological conditions.
Results
We present a simple theoretical framework to show how heterogeneous differentiation in a twomasterregulator paradigm can be governed by a signaling network motif common to all subsets of CD4^{+} T cells. With this motif, a population of naïve CD4^{+} T cells can integrate the signals from their environment to generate a functionally diverse population with robust commitment of individual cells. Notably, two positive feedback loops in this network motif govern three bistable switches, which in turn, give rise to three types of heterogeneous differentiated states, depending upon particular combinations of input signals. We provide three prototype models illustrating how to use this framework to explain experimental observations and make specific testable predictions.
Conclusions
The process in which several types of T helper cells are generated simultaneously to mount complex immune responses upon pathogenic challenges can be highly regulated, and a simple signaling network motif can be responsible for generating all possible types of heterogeneous populations with respect to a pair of master regulators controlling CD4^{+} T cell differentiation. The framework provides a mathematical basis for understanding the decisionmaking mechanisms of CD4^{+} T cells, and it can be helpful for interpreting experimental results. Mathematical models based on the framework make specific testable predictions that may improve our understanding of this differentiation system.
Background
CD4^{+} T helper cells serve as key players in host immune responses by regulating and coordinating a large repertoire of immune cells, such as macrophages, B cells and CD8^{+} T cells. Consequently, CD4^{+} T helper cells are critical in human health ranging from homeostasis to pathogenesis of diseases [1,2]. Central to the functions of CD4^{+} T cells is their ability to produce a wide range of extracellular immunomodulating agents including cytokines and chemokines [3]. In order to correctly direct the immune response to antigen stimulation, CD4^{+} T cells have to secrete appropriate types of cytokines in appropriate amounts, and they achieve this by differentiating into various subtypes of functional CD4^{+} T cells from a pool of precursor cells, known as naïve CD4^{+} T cells. These subsets primarily include T helper 1 (T_{H}1), T helper 2 (T_{H}2), T helper 17 (T_{H}17) and induced regulatory T (iT_{Reg}) cells. Each subtype of CD4^{+} T cells produces a distinctive spectrum of cytokines, and in each of these subtypes there is typically one key transcription factor, or master regulator, that is highly expressed and controls the expression of downstream genes, including those encoding the lineage specific cytokines. The master regulators for the four functional subsets are Tbet, GATA3, RORγt and Foxp3, respectively [3].
The differentiation of CD4^{+} T cells is a highly controlled process, and the lineage specificity of the differentiation process is determined by integrating microenvironmental cues that activate various signaling pathways. These pathways include the T cell receptor (TCR) pathway and the Signal Transducer and Activator of Transcription (STAT) pathways [4,5], which are activated by cognate antigens and cytokines, respectively. Other pathways, such as those associated with Notch and Tolllike receptors (TLRs), are also involved in differentiation of CD4^{+} T cells into distinct lineages [68].
In a few types of chronic infections, the dominance of one subtype of CD4^{+} T cells can be observed [9]. However, most immune responses elicit balanced phenotypes of functional CD4^{+} T cells and their effector molecules, suggesting the importance of maintaining the diversity and flexibility of functional CD4^{+} T cells [10,11]. The importance of balancing the phenotypic composition is further corroborated by the fact that inappropriate dominance of particular subtype(s) of CD4^{+} T cells is often associated with inflammatory disorders [1214]. It is not surprising to observe the balanced phenotypes of CD4^{+} T cells in vivo, given the plausible heterogeneous microenvironments of the naïve CD4^{+} T cells, which may stimulate the differentiation into multiple subtypes of functional CD4^{+} T cells. Interestingly, however, highly purified naïve CD4^{+} T cells can be induced to differentiate into multiple subtypes simultaneously in certain homogeneous in vitro experimental conditions [1521]. Also interesting are the observations that optimum experimental conditions for generating homogeneous subsets of CD4^{+} T cells often include conditions that block the differentiation of undesired subsets [3]. These observations suggest that some highly regulated mechanisms, intrinsic to naïve CD4^{+} T cells, generate and maintain phenotypic heterogeneity of functional CD4^{+} T cells. In vitro assays showing heterogeneous differentiation recapitulate, at least in part, the balanced CD4^{+} T cell populations observed in vivo. Understanding situations of induced heterogeneous differentiation will shed light on the mechanisms controlling the response of populations of CD4^{+} T cells under physiological conditions.
Although the overexpression of one type of master regulator is generally considered the hallmark of the differentiation of one subtype of CD4^{+} T cells, it has been recently discovered that cells highly expressing two types of master regulators exist in vivo[16,17,2226], and some of these 'doublepositive' phenotypes have been shown to be important in responding to pathogens [16,17,26]. Consistent with in vivo studies showing that the numbers of singlepositive and doublepositive CD4^{+} T cells can be increased in comparable proportions upon pathogenic challenges [16], in vitro induction of the differentiation of doublepositive CD4^{+} T cells often requires heterogeneous differentiation, which is accompanied by the differentiation of singlepositive phenotypes [1517]. Some doublepositive CD4^{+} T cells can be generated by reprogramming the singlepositive phenotypes, which also results in a heterogeneous population containing both singlepositive and doublepositive cells [23,24]. These experiments provide us with the clues to the conditions for generating doublepositive phenotypes and highlight the intimate link between the doublepositive phenotype and heterogeneous differentiation.
In most experiments demonstrating induction of heterogeneous differentiation, the expression levels of master regulators controlling two population subsets are examined at the single cell level. Despite the limited scope of these experiments in terms of the number of subsets considered, significant diversity of heterogeneous differentiation has been revealed. In a particular differentiation event, one can obtain one of the following types of heterogeneous populations (Figure 1): a population containing two types of singlepositive cells [18], a population containing one type of singlepositive cells and doublepositive cells [17], and a population containing two types of singlepositive cells and doublepositive cells [15]. The diversity of heterogeneous differentiation in this minimum paradigm might be only the tip of an iceberg of complexity involving heterogeneous differentiation of all subsets of CD4^{+} T cells, but understanding a minimal system with only two classical subtypes is surely the place to start.
Figure 1. Induced heterogeneous differentiation of CD4^{+}T cells with respect to a pair of master regulators (X and Y). a. Diversity of cell phenotypes during induced differentiation. In the undifferentiated cell, the expression level of both X and Y are low. When the cell is differentiated, three possible functional phenotypes can be obtained: X singlepositive cell, Y singlepositive cell and doublepositive cell. b. Three types of induced heterogeneous differentiation. In a differentiation event, a group of naive cells can be differentiated into two types of singlepositive cells (Type 1), one type of singlepositive cell and DP cell (Type 2) or all three functional phenotypes (Type 3).
Previously, mathematical modeling has advanced our understanding of CD4^{+} T cell differentiation [2732]. In particular, Höfer et al. [27] used a mathematical model to explain T_{H}2 cell fate memory created by positive feedbacks in the signaling network; Mariani et al. [28] used a similar model to demonstrate the robust lineage choice between T_{H}1 and T_{H}2 cells; Yates et al. [29] linked the dynamics of master regulators to the phenotypic composition of T_{H}1 and T_{H}2 cells during differentiation and reprogramming; van den Ham et al. [30] used a generic model to describe the switches among all CD4^{+} T cell lineages; and Naldi et al. [32] developed a Booleannetwork model that takes all four lineages of CD4^{+} T cells into consideration. We recently used a mathematical model to study the reciprocal differentiation of T_{H}17 and iT_{Reg} cells, in which heterogeneous differentiation is observed [33]. It is unclear, however, how a broader spectrum of CD4^{+} T cells can be involved in heterogeneous differentiation and what determines the observed types of differentiated states.
Here, we propose a simple theoretical framework for understanding the heterogeneous differentiation of CD4^{+} T cells. We analyze the dynamic properties of a signaling network motif common to all CD4^{+} T cell lineages. We show that, at the level of cell populations, this motif can generate all possible homogeneous and heterogeneous phenotypic compositions with respect to a pair of master regulators, and at the singlecell level it ensures the robust commitment of a particular choice of differentiated state. Two types of positive feedback loops in this network motif govern three types of bistable switches, which in turn, result in three types of heterogeneous differentiation upon receiving appropriate combinations of input signals. This framework facilitates not only an intuitive understanding of the complex process by which CD4^{+} T cells integrate multiple signals to give rise to multiple functional phenotypes, but also the construction of more detailed mathematical models for studying CD4^{+} T cell differentiation. We provide three prototype models illustrating how to use this framework to explain experimental observations and make specific testable predictions.
Results and discussion
A basal signaling network motif is proposed to govern the differentiation of all lineages of CD4^{+} T cells
To consider the heterogeneous differentiation of CD4^{+} T cells, we introduce a minimal model based on a pair of master regulators (proteins X and Y). We neglect the influence of other master regulators during the differentiation process. In the undifferentiated (naïve) cell, the expression levels of X and Y are both low, and the stable expression of either X or Y marks the differentiation event. Three phenotypes can be observed upon differentiation: X singlepositive (XSP) cell, Y singlepositive (YSP) cell, and doublepositive (DP) cell (Figure 1A). In the model, heterogeneous differentiation is defined as the process in which more than one functional (nonnaïve) phenotypes can be observed upon uniform treatment of a population of simulated naïve cells (see Methods).
In this minimum paradigm, three types of heterogeneous differentiation can be induced: 1) two different types of singlepositive cells are differentiated simultaneously from naïve precursors; 2) one type of singlepositive cells differentiates simultaneously with doublepositive cells; and 3) both types of singlepositive cells differentiate simultaneously with double positive cells (Figure 1B). We define these three scenarios as Type 1, 2 and 3 heterogeneous differentiations, respectively.
We next propose a basal network motif that governs cell differentiation in this minimal model. Based on known molecular interactions, we notice that the four master regulators of CD4^{+} T cells are all involved in signaling networks of similar topologies (Figure 2AC). From these examples, we introduce a ‘basal motif’ (Figure 2D). In the basal motif, two master regulators (X and Y) mutually inhibit each other’s expression, while activating their own production. Two types of signals are responsible for activating the expression of the master regulators: a 'primary signal' (S1) which is sufficient to fully upregulate at least one master regulator, and two polarizing signals (S2 and S3) which favor the expression of one master regulator or the other (X and Y, respectively) but are not sufficient to upregulate their expression in the absence of a primary signal (Figure 2D). Each influence relationship in this basal motif has direct biological meaning, but some components in this motif may represent different biological entities in different dualmasterregulator networks. For example, in the T_{H}1T_{H}2 network (Figure 2B) the primary signal represents the TCR ligands, whereas in the iT_{Reg}T_{H}17 network (Figure 2C) it represents a combined treatment of TCR ligands and TGFβ, which is justified by the fact that both TCR and TGFβ signaling pathways activate both Foxp3 and RORγt. Note that the signals, which are treated as parameters in our models, represent exogenous cytokine doses only, not endogenous cytokines produced by T cells upon activation. The latter are represented in part by the autoactivation relations.
Figure 2. Basal network motif controlling heterogeneous differentiation in the two master regulator paradigm. Solid green arrow: activation influence in which the activator alone can switch on the expression of the target protein. Dashed green arrow: activation influence in which the activator alone cannot switch on the expression of the target protein. Red arrow: inhibition influence. Protein name in parenthesis: possible intermediate protein for the positive feedback loop. a. Prototype Model 1: heterogeneous differentiation of T_{H}1 and T_{H}2. b. Prototype Model 2: heterogeneous differentiation of T_{H}1 and T_{H}17. c. Prototype Model 3: heterogeneous differentiation of iT_{Reg} and T_{H}17. d. The basal network motif.
In Table 1, we list the generic signaling components and their corresponding biological entities for each prototype model. Note that a TCR ligand is a typical example of a primary signal, and certain groups of cytokines correspond to polarizing signals. In Table 2, we list the evidences for all molecular influences of each prototype model.
Table 1. Signaling components in basal motif and their corresponding biological components in prototype models
Table 2. Evidences for molecular influences in prototype models
We first analyze Type 1 heterogeneous differentiation using the core motif, in the absence of autoactivation, and then we use the full version of the basal motif to explain all three types of heterogeneous differentiation.
The basal motif without autoactivations can generate Type 1 heterogeneous differentiation
The symmetric case
Consider first the case of perfectly symmetrical parameter settings (Additional file 1: Table S1 Generic Model 1) for the core motif without selfactivations. (See the Methods section for a description of our mathematical model of the signaling motifs.) In the absence of exogenous signals, the system persists in the stable ‘doublenegative’ state corresponding to naïve cells (Figure 3A). Small positive values of the primary signal (0 < S1 < 0.704) drive the expression of modest amounts of both master regulators in a single cell. Larger values (0.704 < S1 < 2.396) destabilize the coexpression state and give rise to two new (alternative) stable steady states: the XhighYlow state and the XlowYhigh state, which correspond to XSP and YSP cells, respectively (Figure 3B). The basins of attraction of these two states are separated by the diagonal line (X = Y) through the state space. When the primary signal is extremely strong (S1 > 2.396), the system is attracted to a unique stable steady state (XhighYhigh), corresponding to a DP cell (Figure 3C). Bifurcation analysis on these steady states shows that the system undergoes pitchfork bifurcations at S1 = 0.704 and at S1 = 2.396 (Figure 3D), a typical type of bifurcation obtained for dynamical systems with perfect symmetry [5254]. Saturation of the primary signal may prevent cells from reaching the DP state (Additional file 2: Figure S1A and B).
Additional file 1. 4 supplementary tables and legends for supplementary figures.
Format: DOC Size: 211KB Download file
This file can be viewed with: Microsoft Word Viewer
Additional file 2. Figure S1. Effects of primary signal saturation.
Format: TIFF Size: 121KB Download file
Figure 3. Analyses of the core motif with symmetrical parameters. ac. Phase plane portraits for three values of primary signal strength (zero, intermediate, high), in the absence of polarizing signals (S2 = S3 = 0). Green curve: X nullcline; red curve: Y nullcline; blue arrow: representative vector in the phase space; closed circle: stable steady state; open circle: unstable steady state; gray curve: separatrix. d. Oneparameter bifurcation diagram for steady state level of X as a function of primary signal S1. Solid curve: stable steady state; dashed curve: unstable steady state. e. Twoparameter bifurcation diagram with respect to primary signal S1 and polarizing signal S2, with S3 = 0. Solid curve: locus of pitchfork bifurcation points. The pitchfork bifurcation points coalesce and disappear at S2 = 0.357. f. Bidirectional twoparameter bifurcation diagram with respect to primary signal S1 and polarizing signals S2 and S3. Top half: S1—S2 diagram, with S3 = 0, as in panel E. Bottom half: S1—S3 diagram, with S2 = 0. The types of stable steady states in each region are annotated as colored circles. Adjoined circles: multistability. See Figure 1 for interpretation of the color scheme. g. Simulation results for treatment of a population of cells simultaneously with primary and polarizing signals. h. Simulation results for sequential treatment: polarizing signal followed by primary signal. I. Simulation results for sequential treatment: primary signal followed by polarizing signal. In GI, the heterogeneity scores with respect to XSP and YSP are plotted.
The presence of a polarizing signal breaks the symmetry of the system, resulting in a pitchfork bifurcation with broken symmetry (Additional file 3: Figure S2A and B). To analyze the influence of polarizing signals on this dynamical system, we plot twoparameter bifurcation diagrams with respect to the primary signal and to each of the polarizing signals (e.g., Figure 3E, for S1 and S2). In Figure 3F we plot a ‘bidirectional’ twoparameter bifurcation diagram, with S2 versus S1 plotted ‘up’ and S3 versus S1 plotted ‘down’ (see Methods for details). In Figure 3F we see a bistable region (bounded by the red curves) for moderate values of the primary signal strength (0.72.3 units) and for low values (0–0.35 units) of either of the polarizing signal strengths. Within the bistable region are found the two types of singlepositive states. Outside the bistable region are found unique steady state solutions that vary continuously from the naïve state on the left to the doublepositive state on the right, through intermediate region (0.7 < S1 < 2.3) dominated by XSP cells (for S2 > 0) or by YSP cells (for S3 > 0). Because of the perfect symmetry of the parameters, both of the cusps of the bistable region lie on the Xaxis.
Additional file 3. Figure S2. Hysteresis effect of the ‘reprogramming’ bistable switch.
Format: TIFF Size: 319KB Download file
In order to predict the response of this regulatory system to changing stimuli (S1 and S2, or S1 and S3), we must be careful in interpreting the effects of trajectories crossing the twoparameter bifurcation diagram in Figure 3F. If we fix the polarizing signals at S3 = 0, S2 = 0.1 and increase the primary signal from 0 to 3, as in Additional file 3: Figure S2A and B, we see that the regulatory system passes smoothly from the naïve state (XlowYlow) to the XSP state (XhighYlow) to the DP state (XhighYhigh). The regulatory system passes over the bistable region without undergoing any abrupt changes of the state (bifurcation) or exhibiting hysteresis effects. On the other hand, if we fix the primary signal at S1 = 1.5 and increase one of the polarizing signals (either S2 or S3), as in Additional file 3: Figure S2 C and D, we see that the regulatory system starts in one of the singlepositive state and jumps abruptly to another singlepositive state at a saddlenode bifurcation point. Also, the system exhibit hysteresis because, if the polarizing signal is reduced to zero after the jump occurs, the regulatory system remains stuck in the stable ‘flipped’ state (XSP if S2 increases/decreases, YSP if S3 increases/decreases). We call this type of response a ‘reprogramming’ switch, because the control system flips irreversibly between alternative singlepositive states. On the contrary, transitions from the naïve or the DP state to either one of the singlepositive states are smooth and reversible (they do not invoke reprogramming).
We next show that this network motif can generate heterogeneous differentiation and identify the parameter region in which a heterogeneous population can be obtained. To this end we simulate the induced differentiation process in a group of cells (with small celltocell variability) exposed to various combinations of primary (S1) and polarizing signals (either S2 or S3). For each combination of S1 and S2 (or S3), we compute the percentages of cells of different phenotypes in the final (steady state) differentiated population. We plot these percentages (as heat maps) over the coordinates of the bidirectional twoparameter bifurcation diagram (see Additional file 4: Figure S3AD). We summarize these results with a ‘heterogeneity score’ (see Methods) to highlight the region of parameter space that can generate heterogeneous populations (Figure 3G). Not surprisingly, in the absence of strong polarizing signals (S2 ≈ 0 and S3 ≈ 0), the primary signal can induce heterogeneous differentiation of two singlepositive phenotypes (Figure 3G, bright area). This is because of the close proximity of the naïve states to the separatrix, and the presence of celltocell variability which can bias individual cells towards different phenotypes (Additional file 4: Figure S3E). The polarizing signal, on the other hand, makes the differentiation into one singlepositive phenotype more likely, which can result in homogeneous differentiation once it is sufficiently strong (Figure 3G, dark area).
Additional file 4. Figure S3. Simulation results for the core motif with symmetrical parameters.
Format: TIFF Size: 731KB Download file
We next explore how the cell population responds to sequential stimuli rather than simultaneous stimuli. If the population is stimulated first by a polarizing signal and then, after the cells have reached their steady states, the simulations are continued in the presence of primary signal, we find that the response to sequential stimuli is very similar to the response to simultaneous stimuli (Figure 3H). But when we switch the sequence of the stimuli, the polarizing signal fails to influence cell fate in the bistable region, resulting in heterogeneous populations in this region (Figure 3I). This is due to a hysteresis effect, which prevents reprogramming by polarizing signals that are insufficiently strong. These results suggest that polarizing signals can influence cell fate determination until the induction of differentiation, after which their influence is greatly reduced.
Broken symmetry
The preceding analysis is based on a set of perfectly symmetrical parameters in the signaling network, although the exogenous polarizing signals can act as ‘symmetry breakers’. How differently does the regulatory system behave if its intrinsic kinetic parameters are not perfectly symmetrical? For illustrative purposes, we use a representative set of asymmetrical parameter values (Additional file 1: Table S1 Generic Model 2). Because of the asymmetries, the primary signal upregulates the two master regulators at different thresholds (Figure 4A and B), and the bistable region of the bidirectional twoparameter bifurcation diagram is reoriented so that its cusps are located on different sides of the Xaxis (Figure 4C). When we stimulate cell populations with combinations of primary and polarizing signals, we find that the parameter region that gives rise to heterogeneous populations is not coincident with the Xaxis. Instead, the ‘heterogeneous’ region forms a patch that intersects the Xaxis (Figure 4D). In this situation, the system requires a specific range of primary signal strength to generate a heterogeneous population. On the other hand, the primary signal now gains some control over cell fate determination, in addition to its ability to trigger the differentiation. For a similar network in B cells, Sciammas et al. [55] recently showed that the strength of the B cell receptor signal (primary signal) can determine cell fate because of the asymmetry of the network.
Figure 4. Analyses of the core motif with asymmetrical parameters. a and b. Oneparameter bifurcation diagram for steady state levels of X and Y as functions of primary signal S1 (S2 = S3 = 0). c. Bidirectional twoparameter bifurcation diagram with respect to primary signal S1 and polarizing signal S2 or S3. See legend of Figure 3 panels D and E for the interpretation of curves and colored circles. df. See legend of Figure 3 Panels GI for simulation conditions.
The effects of sequential stimuli in the asymmetrical model are similar to their effects in the symmetrical model (Figure 4E and F).
Up to this point, we have assumed that the relaxation rates of X and Y are identical . Breaking this symmetry changes the parameter combinations that generate heterogeneous differentiation without changing the bifurcation diagram (Additional file 5: Figure S4). This result, together with the responses to sequential stimuli discussed earlier, shows that although the bistable region is critical to obtaining heterogeneous differentiation, the exact phenotypic composition within the bistable region also depends on the kinetics of the signal inputs and the intrinsic relaxation rates of the master regulators.
Additional file 5. Figure S4. Simulation results with different relaxation rates of X and Y.
Format: TIFF Size: 336KB Download file
We suggest that biological signaling networks of this type (i.e., those resembling the basal motif) may have evolved to take advantage of either symmetrical or asymmetrical types of behavior. A typical asymmetrical design is found in the T_{H}1 and T_{H}2 paradigm, in which TCR signaling not only triggers the heterogeneous differentiation of both T_{H}1 and T_{H}2, but also regulates their phenotypic compositions depending on signal strength (discussed in detail in later section). With this understanding, one can design experiments to study more detailed signalcontrol principles of a particular signaling network governing heterogeneous differentiation.
The basal network motif with additional positive feedback loops can generate all types of heterogeneous differentiation
Previously, mathematical modelers found that interconnected positive feedback loops can give rise to complex multistability in CD4^{+} T cell differentiation [28] and elsewhere [54]. It is still not clear, however, how these different multistable regions depend on the interconnection of multiple positive feedback loops, nor how one can use biologically relevant signals to guide cells into various multistable regions, where heterogeneous differentiation might occur. In this section, we show that our basal motif can give rise to complex multistability, we clarify the effects of the additional positive feedback loops using bifurcation analysis, and we explain the biological meaning of each parameter region in the context of the heterogeneous differentiation of CD4^{+} T cells.
For illustrative purpose, we first choose another set of perfectly symmetrical parameters (Additional file 1: Table S1 Generic Model 3). This model differs from Generic Model 1 in that the doublenegative feedback (mutual inhibition) is not strong enough to create bistability. Nonetheless, with the addition of symmetrical increase of autoactivation loops, a bistable region first appears in the intermediate range (1.7 < S1 < 2.4) of the primary signal (Additional file 6: Figure S5A), similar to the case of Generic Model 1 (Figure 3D). Further increase of the autoactivation weights enlarges the bistable region, and at a critical point (weights = 1.8), the pitchfork bifurcation changes from supercritical (Additional file 6: Figure S5A, weights = 1.5) to subcritical (Additional file 6: Figure S5B, weights = 3.2). Beyond the transition from supercritical to subcritical, each pitchfork bifurcation gives rise to two saddlenode bifurcation points (Additional file 6: Figure S5B and C). On the bidirectional (S1S2S3) twoparameter bifurcation diagram (Figure 5A), each cusp region 'folds back' to form three interconnected cusp regions, which govern two new bistable regions and one tristable region (Figure 5A). Further increase of the autoactivation weights enlarges the original bistable region as well as the newly formed multistable regions. Eventually, the plane on the bidirectional twoparameter bifurcation diagram is divided into 11 regions with distinct stability features (Figure 5B).
Additional file 6. Figure S5. Additional bifurcation analyses of the full basal motif.
Format: TIFF Size: 508KB Download file
Figure 5. Analyses of the basal motif with autoactivation relations. a. Bidirectional twoparameter bifurcation diagram with respect to primary signal S1 and polarizing signals S2 and S3 for intermediate weight of autoactivation relations (ω = 1.8). Insets show the zoomedin view of the cusp regions of the bistable region. b. Bidirectional twoparameter bifurcation diagram with respect to primary signal S1 and polarizing signals S2 and S3 for high weight of autoactivation relations (ω = 3.2). The types of stable steady states in each region are annotated as colored circles. Adjoined circles: multistability. See Figure 1 for interpretation of the color scheme. Light blue area: bistable region governing differentiation switch. Light green area: bistable region governing reprogramming switch. Light yellow area: bistable region governing coexpression switch. cf. Various types of heterogeneity scores are plotted for high weight of autoactivation relations (ω = 3.2). c. The heterogeneity scores with respect to XSP and YSP. d. The heterogeneity scores with respect to XSP and DP. e. The heterogeneity scores with respect to YSP and DP. f. The heterogeneity scores with respect to XSP, YSP and DP.
We clarify this unique twoparameter bifurcation diagram as follows. If the autoactivation loops are absent or weaker, the parameter region outside of the reprogramming switch bistable region (Figure 3F) is continuous and monostable, although it can represent four types of steady states. Essentially, strong autoactivation loops create folding in this monostable region so that it is divided into four monostable regions separated by four new bistable regions. This structure effectively creates an additional level of robustness of cell fate commitment, which is rendered by two new types of bistable switches, in addition to the reprogramming switch. One type of switch consists of the two bistable regions located at lower range of the primary signal (Figure 5B, light blue areas), which controls differentiation/dedifferentiation commitment, i.e. the switches from or to the naïve state (Additional file 6: Figure S5D and E). Another type of switch consists of the two bistable regions located at higher range of the primary signal (Figure 5B, light yellow areas), which controls coexpression commitment, i.e. the switches from or to the doublepositive state (Additional file 6: Figure S5D and E). We define these two switches as the ‘differentiation switch’ and the ‘coexpression switch’ respectively. The tristable regions in this diagram are the overlapping areas between the bistable regions governed by the reprogramming switch and either the differentiation or the coexpression switch. In fact, extremely high weights (>4) of autoactivation may give rise to a tetrastable region, where the three types of the bistable regions overlap (Additional file 6: Figure S5C).
In summary, the positive feedback loop involving mutual inhibition of the master regulators can create the reprogramming switch, and additional feedback loops involving autoactivation can enhance the robustness of the reprogramming switch and create the differentiation switch and the coexpression switch. The features of the three bistable switches are listed in Table 3.
Table 3. Features of three bistable switches obtained with the basal motif
We next ran simulations to check whether these regions of multistability are correlated to various types of heterogeneous differentiation. Our results show that Type 1 heterogeneous differentiation can be induced in the reprogramming switch region (Figure 5C) (this is consistent with the results obtained with the core motif), Type 2 heterogeneous differentiation can be induced in the coexpression bistable switch regions (Figure 5D and E), and Type 3 heterogeneous differentiation can be induced in the tristable region consisting of three functional (nonnaïve) states (Figure 5F). These types of heterogeneous differentiations are all robust in terms of single cell commitment because the corresponding parameter regions admit a variety of stable steady states.
Positive feedback loops have long been recognized as mechanisms for biological switches [5658]. We have demonstrated that two types of positive feedback in the CD4^{+} T cell differentiation network underlie three types of bistable switches that govern the transitions among different phenotypes of those T cells. In addition to ensuring the robust commitment, the multistability created by positive feedback loops may be used to generate phenotypic diversities of various types. In this context, the biological functions of the positive feedback loops are seen as more versatile than giving rise to simple onoroff switches.
Our theoretical analysis of the basal regulatory motif (Figure 2D) started with symmetrical parameter values and then considered the effects of broken symmetries. In the next section, we show how nonsymmetrical prototype models of heterogeneous differentiation among real lines of CD4^{+} T cells can be studied within this unifying framework despite their diverse features.
Mathematical models based on the theoretical framework can be used to understand experimental results and make testable predictions
In this section we discuss three prototype models for studying heterogeneous differentiation of CD4^{+} T cells. The first two models are aimed to explain some interesting biological phenomena that were not studied previously with mathematical modeling. The third one is a simplified version of our previous model [33], but we have made it more accessible by using the framework presented here. Because of their limited scope, none of these models are intended to provide a comprehensive understanding of the corresponding biological systems. Rather, our intention is to illustrate how to use the modeling framework to explain observed heterogeneous differentiation and make testable predictions.
Prototype Model 1: Heterogeneous differentiation of T_{H}1 and T_{H}2 cells
Previous mathematical models successfully described the dynamic behavior and the underlying molecular control system of the reciprocal differentiation of T_{H}1 and T_{H}2 cells [2731]. However, heterogeneous differentiation of T_{H}1 and T_{H}2 cells and its underlying molecular controls were not studied with these models. Yamashita et al. [18] discovered that the heterogeneous differentiation of T_{H}1 and T_{H}2 cells can be obtained with antigenic stimulations. Similar observations were obtained by Hosken et al. [20], and Messi et al. [21]. We have built a mathematical model, based on the influence diagram in Figure 2A, to describe heterogeneous differentiation of T_{H}1 and T_{H}2 cells. The parameter values for the model are listed in Additional file 1: Table S2.
Figure 6A shows the bidirectional twoparameter bifurcation diagram, and Figure 6B shows the simulation results as the heterogeneity score with respect to the two singlepositive phenotypes. Our simulation results suggest that exogenous polarizing signals, i.e. IL4 and IL12, are not sufficient to trigger differentiation. They must be accompanied by a sufficiently high dose of antigenic stimulant (TCR signal) to trigger the differentiation into the corresponding phenotypes. This conclusion is in agreement with previous experimental results [18]. High strength of TCR signal alone (>1 unit) or with intermediate level of IL4 (0.3 unit) was sufficient to induce the differentiation of two singlepositive phenotypes. With increasing strengths of TCR signal, our simulations show a spectrum of heterogeneous populations with increasing percentages of T_{H}2 cells and decreasing percentage of T_{H}1 cells. The following experimental findings are consistent with our simulation. Messi et al. [21] observed the heterogeneous differentiation of T_{H}1 and T_{H}2 with IL4 and antigenic stimulant. Yamashita et al. [18] observed a similar pattern of heterogeneous populations with increasing doses of antigenic stimulant in the presence of an intermediate level of IL4. Hosken et al. [20] also observed such pattern with a different type of antigenic stimulant, although only a narrow range of stimulant concentrations could give rise to heterogeneous populations. Clearly, our model predicts that in order to achieve comparable proportions of T_{H}1 cells and T_{H}2 cells, one would need a higher dose of antigenic stimulant without exogenous IL4 as compared to with exogenous IL4. Based on the bifurcation diagram, we also predict that a slow increase of stimulant concentration favors the differentiation of T_{H}1 cells. Additionally, the simulation results and bifurcation analysis show that the doublepositive phenotype can be obtained in the presence of T_{H}1 polarizing signals. Hegazy et al. [24] have discovered that exogenous T_{H}1 polarizing signals can reprogram T_{H}2 cells into Tbet^{+}GATA3^{+} cells in the presence of antigenic stimulant. Our model predicts that the differentiation of such doublepositive phenotype can be directly induced by high dose of antigenic stimulant (>2 units) in the presence of exogenous T_{H}1 polarizing signals (0.5 unit), and the differentiation is likely to be heterogeneous with the concurrent induction of two types of singlepositive cells, in addition to the doublepositive cells. If we reduce the autoactivation weight of GATA3 (see Methods), then the TCR signal primarily triggers the differentiation of T_{H}1 cells instead of a heterogeneous population (Figure 6C and D). Maruyama et al. [59] demonstrated that TCR signal alone can induce a significant fraction of GATA3^{+} cells (this is consistent with the experimental findings mentioned above), and blocking the autoactivation feedback between GATA3 and IL4 prevents the induction of GATA3^{+} cells. Our model predicts that the population may be dominated by T_{H}1 cells under this condition.
Figure 6. Analyses of Prototype Model 1 (heterogeneous differentiation of T_{H}1 and T_{H}2 cells). a. Bidirectional twoparameter bifurcation diagram with respect to primary signal TCR and polarizing signals IL12 and IL4. b. Simulation results for induced differentiation. The heterogeneity scores with respect to Tbet singlepositive phenotype and GATA3 singlepositive phenotype are shown. c. Same legend as Panel A. The GATA3 autoactivation relation is blocked in the model. d. Same legend as Panel B. The GATA3 autoactivation relation is blocked in the model. In Panels A and C: Adjoined circles: multistability. Blue circle: naïve phenotype. Green circle: Tbet singlepositive phenotype. Red circle: GATA3 single positive phenotype. Yellow: DP phenotype.
Table 4 summarizes the published observations consistent with our simulation results and new predictions based on the bifurcation analyses and simulation results.
Table 4. Summary of simulation results of Prototype Model 1
Prototype Model 2: Heterogeneous differentiation of T_{H}1 and T_{H}17 cells
We build a prototype model to study the heterogeneous differentiation of T_{H}1 and T_{H}17 cells that was recently demonstrated by Ghoreschi et al. [17]. The influence diagram of the model is shown in Figure 2B, and the parameter values are listed in Additional file 1: Table S3. In the presence of TCR signal alone, the simulated population is dominated by T_{H}1 cells (Figure 7A and B). When the TCR signal is combined with IL23 + IL1 polarizing signal, the induced population contains both the Tbet^{+}RORγt^{} singlepositive phenotype and the Tbet^{+}RORγt^{+} double positive phenotype (Figure 7A and B). When the TCR signal is combined with TGFβ (another polarizing signal), the population is dominated by the Tbet^{}RORγt^{+} singlepositive phenotype (Figure 7C and D). These results are consistent with the observations of Ghoreschi et al. [17]. Our model predicts that lowering the TCR signal strength may result in the reprogramming from Tbet^{+}RORγt^{+} double positive phenotype to Tbet^{+}RORγt^{} single positive phenotype even in the presence of a strong IL23 + IL1 signal and that when low dose of TGFβ + IL6 (≈0.4 unit) is used, one may observe the heterogeneous differentiation of T_{H}1 and T_{H}17 cells. Also, the model recapitulates the scenario in which knocking out Tbet genes resulted in the homogeneous differentiation into Tbet^{}RORγt^{+} singlepositive phenotype when either of the polarizing signals is used (Additional file 7: Figure S6) [17].
Additional file 7. Figure S6. Simulation results of Prototype Model 2 (heterogeneous differentiation of T_{H}1 and T_{H}17 cells) with Tbet knockedout.
Format: TIFF Size: 596KB Download file
Figure 7. Analyses of Prototype Model 2 (heterogeneous differentiation of T_{H}1 and T_{H}17 cells). a. Twoparameter bifurcation diagram with respect to primary signal TCR and polarizing signal IL23 + IL1. b. Simulation results for induced differentiation. The heterogeneity scores with respect to Tbet singlepositive phenotype and DP phenotype are shown. c. Twoparameter bifurcation diagram with respect to primary signal TCR and polarizing signal TGFβ + IL6. d. Simulation results for induced differentiation. Heterogeneity scores with respect to Tbet singlepositive phenotype and RORγt singlepositive phenotype are shown. In Panels A and C: Adjoined circles: multistability. Blue circle: naïve phenotype. Green circle: Tbet singlepositive phenotype. Red circle: RORγt singlepositive phenotype. Yellow: DP phenotype.
Simulation results with testable predictions are summarized in Table 5.
Table 5. Summary of simulation results of Prototype Model 2
Prototype Model 3: Heterogeneous differentiation of iT_{Reg} and T_{H}17 cells
Heterogeneous differentiation of iT_{Reg} and T_{H}17 cells has been observed in many experiments [15,16,19]. Here we present a prototype model based on the influence diagram (Figure 2C) and the parameter values (Additional file 1: Table S4). The model shows that a combination of TGFβ and TCR signal can drive a heterogeneous population containing Foxp3^{+}RORγt^{}, Foxp3^{}RORγt^{+} and Foxp3^{+}RORγt^{+} phenotypes (Figure 8A and B, tristable region at TCR + TGFβ signal ≈ 1.8). Raising the strength of TGFβ + TCR signal or adding IL6 (a T_{H}17 polarizing signal) can skew the population into Foxp3^{}RORγt^{+} and Foxp3^{+}RORγt^{+} phenotypes (Figure 8A and B, bistable region in the upper plot at highest level of TCR + TGFβ signal). These results are in agreement with previous experimental observations [15,16]. Predictions made from the model include: 1) an intermediate TGFβ + TCR signal (1–1.5 units) favors heterogeneous differentiation of Foxp3^{+}RORγt^{} and Foxp3^{}RORγt^{+} populations; 2) an intermediate level of TGFβ + TCR signal (1–1.5 units) with an iT_{Reg} polarizing signal produces a homogeneous Foxp3^{+}RORγt^{} population; and 3) a high level of TGFβ + TCR signal (>2 units) with an iT_{Reg} polarizing signal induces heterogeneous Foxp3^{+}RORγt^{} and Foxp3^{+}RORγt^{+} populations.
Figure 8. Analyses of Prototype Model 3 (heterogeneous differentiation of iT_{Reg}and T_{H}17 cells). a. Bidirectional twoparameter bifurcation diagram with respect to primary signal TCR + TGFβ and polarizing signals ATRA/IL2 and IL6. Adjoined circles: multistability. Blue circle: naïve phenotype. Green circle: Foxp3 singlepositive phenotype. Red circle: RORγt singlepositive phenotype. Yellow: DP phenotype. b. Simulation results of induced differentiation. The heterogeneity scores with respect to Foxp3 singlepositive phenotype, RORγt singlepositive phenotype and DP phenotype are shown.
Simulation results with testable predictions are summarized in Table 6.
Table 6. Summary of simulation results of Prototype Model 3
Conclusions
In this study, we have demonstrated that a simple signaling network motif can be responsible for generating all possible types of heterogeneous populations with respect to a pair of master regulators controlling CD4^{+} T cell differentiation. We showed how naïve CD4^{+} T cells can integrate multiple types of signals to differentiate into populations of diverse phenotypes. We illustrate the theoretical framework with three specific cases and made testable predictions.
It is becoming evident that certain signals can drive the differentiation of multiple lineages of T cells, whereas other environmental cues can skew the outcome to specific phenotypes [60]. Because the proposed basal motif appears commonly in the signaling networks controlling CD4^{+} T cell differentiation, biological examples of this framework are clearly not limited to the prototype models we presented here. For example, it has been recently demonstrated that STAT3 activation is required for T_{H}2 differentiation [61]. This gives the possibility that IL6, which upregulates RORγt via STAT3 activation [62], can act as a primary signal giving rise to heterogeneous T_{H}2 and T_{H}17 populations if the cells are primed with certain amount of other signals, such as TCR, TGFβ and IL4.
Our study suggests the importance of regulated celltocell variations that can be exploited to generate phenotypic diversity in CD4^{+} T cells. The significance of such variations in some other biological systems has been highlighted by other groups. Feinerman et al. [63] discovered that the celltocell variations in the expression levels of some key coreceptors in CD8^{+} T cells can be critical for achieving diversity in TCR responses. Similarly, Chang et al. [64] demonstrated that variations in the expression of stem cell markers can influence the fate of the cell. We have used a simple generic form to account for celltocell variability in this study (i.e. parametric variations), it would be interesting to study which specific variable factors in naïve CD4^{+} T cells can be predictive of the phenotypic compositions in an induced population. Harnessing such factors might be useful for finetuning the immune system to prevent and treat diseases.
Our modeling approach has the advantage of describing nonlinear responses in biochemical reactions without knowing detailed biochemical mechanisms and kinetics, which are generally unavailable for T cell differentiation. It has the disadvantage that parameters in the equations are phenomenological and cannot be related to biochemical reaction rate constants. We expect that other modeling approaches, such as ordinary differential equations with Hill function nonlinearities, will produce results similar to ours.
We are aware of the following limitations of this framework. First, all master regulators of CD4^{+} T cell may influence each other during differentiation. Thus considering only a pair of master regulators may not be sufficient to describe all important components governing the heterogeneous differentiation of CD4^{+} T cells. Secondly, celltocell communication is neglected in our models of cell population. We assume that our models describe the initial phase of differentiation and that the phenotypic compositions of the population do not change significantly during the differentiation process. The validity of this assumption needs to be examined in future studies.
Methods
Dynamical model
We modeled the signaling network motifs with a generic form of ordinary differential equations (ODEs) that describe both gene expression and protein interaction networks [6567]. Each ODE in our model has the form:
Where X_{i} is the activity or concentration of protein i. On a time scale = 1/γ_{i}, X_{i}(t) relaxes toward a value determined by the sigmoidal function, F, which has a steepness set by . The basal value of F, in the absence of any influencing factors, is determined by . The coefficients determine the influence of protein j on protein i. N is the total number of proteins in the network.
All variables and parameters are dimensionless. One time unit in our simulations corresponds to 1.5 days. Parameter values are listed in supplementary tables.
All simulations and bifurcation analyses were performed with PyDSTool, a software environment for dynamical systems [68].
Bifurcation diagrams
In order to visualize the response of the T cell differentiation network to multiple signals (a primary differentiation signal and two types of polarizing signals), we have employed bidirectional twoparameter bifurcation diagrams, as in [69]. The two twoparameter bifurcation diagrams share the same primary bifurcation parameter (the primary differentiation signal, S1) on the horizontal axis. The secondary bifurcation parameters (the polarizing signals, S2 and S3) are plotted on the vertical axis: one in the upward direction and the other in the downward direction. The bidirectional twoparameter bifurcation diagram allows one to analyze the response of the regulatory system to the primary signal alone or in combination with either of the polarizing signals. Although this twodimensional representation does not allow a full analysis of the responses to all three types of signals simultaneously, it is very useful in understanding the complex interplay between signals and responses in these heterogeneous differentiation systems. We ran simulations for a population of naïve CD4^{+} T cells, and we overlaid the simulation results on the bidirectional twoparameter bifurcation diagrams, allowing one to visualize the bifurcation analyses and simulation results simultaneously (detailed below).
Celltocell variability
To account for celltocell 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 allowed all of the parameters in our model to change simultaneously, and 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. In our bifurcation analysis of the dynamical system, we considered an imaginary cell that adopts the basal value for each of its parameters, and we defined 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.
In order to simulate the induced differentiation process, we first solved the ODEs numerically with some small initial values of master regulator concentrations in the absence of any exogenous signals. After a short period of time, each simulated cell will find its own, stable ‘doublenegative’ steady state, corresponding to a naïve CD4^{+} T cell. Next, we changed the primary and/or polarizing signals to certain positive values and continued the numerical simulation. If needed, we continued the simulation again with a second change of primary and/or polarizing signals. 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. We repeated this simulation 200 times for a given set of exogenous signals to represent the responses of 200 cells in a population. We made the simple definition that a protein is expressed when its level is greater than 0.5 units. The simulations for a cell population were repeated 40x40 times with primary and polarizing signals of various strengths, and we overlaid the final steady state phenotypic composition on the point with corresponding coordinates on the bidirectional twoparameter bifurcation diagram.
Mutant simulation
The experiment of knocking out GATA3IL4 feedback was simulated with reduced weight of autoactivation of GATA3 to onetenth of the original value. The experiment of knocking out Tbet genes was simulated by setting = −17 (10 times its value in the basal model).
Heterogeneity score
To summarize simulations results with multiple phenotypes and to highlight heterogeneous and homogeneous populations in parameter space, we compute a ‘heterogeneity score’ for a simulation as follows.
The scoring function takes a list of ‘phenotypes of interest’ (), and computes the sum of the pairwise heterogeneities, which are based on the numbers of cells of any two different phenotypes ( and ). The score is normalized with respect to the number of phenotypes of interest (n) and the total number of cells in the population (N). S_{H} ≈ 1 when there are comparable numbers of cells of the phenotypes of interest in the population, S_{H} ≈ −1 when the population is dominated by one phenotype out of all the phenotypes of interest, and S_{H} ≈ 0 when there are few cells with the phenotypes of interest in the population, or the degree of heterogeneity is moderate.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
Conceived and designed the experiments: TH JX LL JJT. Performed the experiments: TH. Analyzed the data: TH JX LL JJT. Wrote the paper: TH LL JJT. All authors read and approved the final manuscript.
Acknowledgements
TH was supported by a Transdisciplinary Team Science Fellowship from the Virginia Bioinformatics Institute, and by an NIH Grant (R01 GM07898905) to JJT. LL was partially supported by grants from NIH (R0164414) and the American Heart Association. JX was supported by an NSF grant (DMS0969417) and an NIAID grant 1R03AI099120. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
References

Dranoff G, Jaffee E, Lazenby A, Golumbek P, Levitsky H, Brose K, Jackson V, Hamada H, Pardoll D, Mulligan RC: Vaccination with irradiated tumor cells engineered to secrete murine granulocytemacrophage colonystimulating factor stimulates potent, specific, and longlasting antitumor immunity.
Proc Natl Acad Sci U S A 1993, 90:35393543. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Moir S, Chun TW, Fauci AS: Pathogenic mechanisms of HIV disease.
Annu Rev Pathol 2011, 6:223248. PubMed Abstract  Publisher Full Text

Zhu J, Yamane H, Paul WE: Differentiation of effector CD4 T cell populations.
Annu Rev Immunol 2010, 28:445489. PubMed Abstract  Publisher Full Text

Nakayama T, Yamashita M: The TCRmediated signaling pathways that control the direction of helper T cell differentiation.
Semin Immunol 2010, 22:303309. PubMed Abstract  Publisher Full Text

O'Shea JJ, Lahesmaa R, Vahedi G, Laurence A, Kanno Y: Genomic views of STAT function in CD4+ T helper cell differentiation.
Nat Rev Immunol 2011, 11:239250. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fang TC, YashiroOhtani Y, Del Bianco C, Knoblock DM, Blacklow SC, Pear WS: Notch directly regulates Gata3 expression during T helper 2 cell differentiation.
Immunity 2007, 27:100110. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Reynolds JM, Pappu BP, Peng J, Martinez GJ, Zhang Y, Chung Y, Ma L, Yang XO, Nurieva RI, Tian Q, Dong C: Tolllike receptor 2 signaling in CD4(+) T lymphocytes promotes T helper 17 responses and regulates the pathogenesis of autoimmune disease.
Immunity 2010, 32:692702. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Amsen D, Blander JM, Lee GR, Tanigaki K, Honjo T, Flavell RA: Instruction of distinct CD4 T helper cell fates by different notch ligands on antigenpresenting cells.
Cell 2004, 117:515526. PubMed Abstract  Publisher Full Text

Sher A, Coffman RL: Regulation of immunity to parasites by T cells and T cellderived cytokines.
Annu Rev Immunol 1992, 10:385409. PubMed Abstract  Publisher Full Text

O'Shea JJ, Paul WE: Mechanisms underlying lineage commitment and plasticity of helper CD4+ T cells.
Science 2010, 327:10981102. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Murphy KM, Stockinger B: Effector T cell plasticity: flexibility in the face of changing circumstances.
Nat Immunol 2010, 11:674680. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Littman DR, Rudensky AY: Th17 and regulatory T cells in mediating and restraining inflammation.
Cell 2010, 140:845858. PubMed Abstract  Publisher Full Text

Romagnani S: Immunologic influences on allergy and the TH1/TH2 balance.
J Allergy Clin Immunol 2004, 113:395400. PubMed Abstract  Publisher Full Text

Mauri C, Williams RO, Walmsley M, Feldmann M: Relationship between Th1/Th2 cytokine patterns and the arthritogenic response in collageninduced arthritis.
Eur J Immunol 1996, 26:15111518. PubMed Abstract  Publisher Full Text

Zhou L, Lopes JE, Chong MM, Ivanov II, Min R, Victora GD, Shen Y, Du J, Rubtsov YP, Rudensky AY, et al.: TGFbetainduced Foxp3 inhibits T(H)17 cell differentiation by antagonizing RORgammat function.
Nature 2008, 453:236240. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lochner M, Peduto L, Cherrier M, Sawa S, Langa F, Varona R, Riethmacher D, SiTahar M, Di Santo JP, Eberl G: In vivo equilibrium of proinflammatory IL17+ and regulatory IL10+ Foxp3+ RORgamma t + T cells.
J Exp Med 2008, 205:13811393. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Ghoreschi K, Laurence A, Yang XP, Tato CM, McGeachy MJ, Konkel JE, Ramos HL, Wei L, Davidson TS, Bouladoux N, et al.: Generation of pathogenic T(H)17 cells in the absence of TGFbeta signalling.
Nature 2010, 467:967971. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yamashita M, Kimura M, Kubo M, Shimizu C, Tada T, Perlmutter RM, Nakayama T: T cell antigen receptormediated activation of the Ras/mitogenactivated protein kinase pathway controls interleukin 4 receptor function and type2 helper T cell differentiation.
Proc Natl Acad Sci U S A 1999, 96:10241029. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Molinero LL, Miller ML, Evaristo C, Alegre ML: High TCR stimuli prevent induced regulatory T cell differentiation in a NFkappaBdependent manner.
J Immunol 2011, 186:46094617. PubMed Abstract  Publisher Full Text

Hosken NA, Shibuya K, Heath AW, Murphy KM, O'Garra A: The effect of antigen dose on CD4+ T helper cell phenotype development in a T cell receptoralpha betatransgenic model.
J Exp Med 1995, 182:15791584. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Messi M, Giacchetto I, Nagata K, Lanzavecchia A, Natoli G, Sallusto F: Memory and flexibility of cytokine gene expression as separable properties of human T(H)1 and T(H)2 lymphocytes.
Nat Immunol 2003, 4:7886. PubMed Abstract  Publisher Full Text

Ayyoub M, Deknuydt F, Raimbaud I, Dousset C, Leveque L, Bioley G, Valmori D: Human memory FOXP3+ Tregs secrete IL17 ex vivo and constitutively express the T(H)17 lineagespecific transcription factor RORgamma t.
Proc Natl Acad Sci U S A 2009, 106:86358640. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Voo KS, Wang YH, Santori FR, Boggiano C, Arima K, Bover L, Hanabuchi S, Khalili J, Marinova E, Zheng B, et al.: Identification of IL17producing FOXP3+ regulatory T cells in humans.
Proc Natl Acad Sci U S A 2009, 106:47934798. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hegazy AN, Peine M, Helmstetter C, Panse I, Frohlich A, Bergthaler A, Flatz L, Pinschewer DD, Radbruch A, Lohning M: Interferons direct Th2 cell reprogramming to generate a stable GATA3(+)Tbet(+) cell subset with combined Th2 and Th1 cell functions.
Immunity 2010, 32:116128. PubMed Abstract  Publisher Full Text

Zhu J, Paul WE: Peripheral CD4+ Tcell differentiation regulated by networks of cytokines and transcription factors.
Immunol Rev 2010, 238:247262. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

AbromsonLeeman S, Bronson RT, Dorf ME: Encephalitogenic T cells that stably express both Tbet and ROR gamma t consistently produce IFNgamma but have a spectrum of IL17 profiles.
J Neuroimmunol 2009, 215:1024. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hofer T, Nathansen H, Lohning M, Radbruch A, Heinrich R: GATA3 transcriptional imprinting in Th2 lymphocytes: a mathematical model.
Proc Natl Acad Sci U S A 2002, 99:93649368. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mariani L, Lohning M, Radbruch A, Hofer T: Transcriptional control networks of cell differentiation: insights from helper T lymphocytes.
Prog Biophys Mol Biol 2004, 86:4576. PubMed Abstract  Publisher Full Text

Yates A, Callard R, Stark J: Combining cytokine signalling with Tbet and GATA3 regulation in Th1 and Th2 differentiation: a model for cellular decisionmaking.
J Theor Biol 2004, 231:181196. PubMed Abstract  Publisher Full Text

van den Ham HJ, de Boer RJ: From the twodimensional Th1 and Th2 phenotypes to highdimensional models for gene regulation.
Int Immunol 2008, 20:12691277. PubMed Abstract  Publisher Full Text

Mendoza L, Xenarios I: A method for the generation of standardized qualitative dynamical systems of regulatory networks.
Theor Biol Med Model 2006, 3:13. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Naldi A, Carneiro J, Chaouiya C, Thieffry D: Diversity and plasticity of Th cell types predicted from regulatory network modelling.
PLoS Comput Biol 2010, 6:e1000912. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hong T, Xing J, Li L, Tyson JJ: A mathematical model for the reciprocal differentiation of T helper 17 cells and induced regulatory T cells.
PLoS Comput Biol 2011, 7:e1002122. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Placek K, Gasparian S, Coffre M, Maiella S, Sechet E, Bianchi E, Rogge L: Integration of distinct intracellular signaling pathways at distal regulatory elements directs Tbet expression in human CD4+ T cells.
J Immunol 2009, 183:77437751. PubMed Abstract  Publisher Full Text

Yamane H, Zhu J, Paul WE: Independent roles for IL2 and GATA3 in stimulating naive CD4+ T cells to generate a Th2inducing cytokine environment.
J Exp Med 2005, 202:793804. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zhu J, Guo L, Watson CJ, HuLi J, Paul WE: Stat6 is necessary and sufficient for IL4's role in Th2 differentiation and cell expansion.
J Immunol 2001, 166:72767281. PubMed Abstract  Publisher Full Text

Usui T, Preiss JC, Kanno Y, Yao ZJ, Bream JH, O'Shea JJ, Strober W: Tbet regulates Th1 responses through essential effects on GATA3 function rather than on IFNG gene acetylation and transcription.
J Exp Med 2006, 203:755766. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zhu J, Jankovic D, Grinberg A, Guo L, Paul WE: Gfi1 plays an important role in IL2mediated Th2 cell expansion.
Proc Natl Acad Sci U S A 2006, 103:1821418219. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mullen AC, High FA, Hutchins AS, Lee HW, Villarino AV, Livingston DM, Kung AL, Cereb N, Yao TP, Yang SY, Reiner SL: Role of Tbet in commitment of TH1 cells before IL12dependent selection.
Science 2001, 292:19071910. PubMed Abstract  Publisher Full Text

Ouyang W, Lohning M, Gao Z, Assenmacher M, Ranganath S, Radbruch A, Murphy KM: Stat6independent GATA3 autoactivation directs IL4independent Th2 development and commitment.
Immunity 2000, 12:2737. PubMed Abstract  Publisher Full Text

Yang XO, Nurieva R, Martinez GJ, Kang HS, Chung Y, Pappu BP, Shah B, Chang SH, Schluns KS, Watowich SS, et al.: Molecular antagonism and plasticity of regulatory and inflammatory T cell programs.
Immunity 2008, 29:4456. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yang L, Anderson DE, BaecherAllan C, Hastings WD, Bettelli E, Oukka M, Kuchroo VK, Hafler DA: IL21 and TGFbeta are required for differentiation of human T(H)17 cells.
Nature 2008, 454:350352. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gorelik L, Constant S, Flavell RA: Mechanism of transforming growth factor betainduced inhibition of T helper type 1 differentiation.
J Exp Med 2002, 195:14991505. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lazarevic V, Chen X, Shim JH, Hwang ES, Jang E, Bolm AN, Oukka M, Kuchroo VK, Glimcher LH: Tbet represses T(H)17 differentiation by preventing Runx1mediated activation of the gene encoding RORgammat.
Nat Immunol 2011, 12:96104. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mukasa R, Balasubramani A, Lee YK, Whitley SK, Weaver BT, Shibata Y, Crawford GE, Hatton RD, Weaver CT: Epigenetic instability of cytokine and transcription factor gene loci underlies plasticity of the T helper 17 cell lineage.
Immunity 2010, 32:616627. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gutcher I, Donkor MK, Ma Q, Rudensky AY, Flavell RA, Li MO: Autocrine transforming growth factorbeta1 promotes in vivo Th17 cell differentiation.
Immunity 2011, 34:396408. PubMed Abstract  Publisher Full Text

Kimura A, Kishimoto T: IL6: regulator of Treg/Th17 balance.
Eur J Immunol 2010, 40:18301835. PubMed Abstract  Publisher Full Text

Elias KM, Laurence A, Davidson TS, Stephens G, Kanno Y, Shevach EM, O'Shea JJ: Retinoic acid inhibits Th17 polarization and enhances FoxP3 expression through a Stat3/Stat5 independent signaling pathway.
Blood 2008, 111:10131020. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mucida D, Park Y, Kim G, Turovskaya O, Scott I, Kronenberg M, Cheroutre H: Reciprocal TH17 and regulatory T cell differentiation mediated by retinoic acid.
Science 2007, 317:256260. PubMed Abstract  Publisher Full Text

Gavin MA, Rasmussen JP, Fontenot JD, Vasta V, Manganiello VC, Beavo JA, Rudensky AY: Foxp3dependent programme of regulatory Tcell differentiation.
Nature 2007, 445:771775. PubMed Abstract  Publisher Full Text

Burgler S, Mantel PY, Bassin C, Ouaked N, Akdis CA, SchmidtWeber CB: RORC2 is involved in T cell polarization through interaction with the FOXP3 promoter.
J Immunol 2010, 184:61616169. PubMed Abstract  Publisher Full Text

Laslo P, Spooner CJ, Warmflash A, Lancki DW, Lee HJ, Sciammas R, Gantner BN, Dinner AR, Singh H: Multilineage transcriptional priming and determination of alternate hematopoietic cell fates.
Cell 2006, 126:755766. PubMed Abstract  Publisher Full Text

Huang S, Guo YP, May G, Enver T: Bifurcation dynamics in lineagecommitment in bipotent progenitor cells.
Dev Biol 2007, 305:695713. PubMed Abstract  Publisher Full Text

Guantes R, Poyatos JF: Multistable decision switches for flexible control of epigenetic differentiation.
PLoS Comput Biol 2008, 4:e1000235. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Sciammas R, Li Y, Warmflash A, Song Y, Dinner AR, Singh H: An incoherent regulatory network architecture that orchestrates B cell diversification in response to antigen signaling.
Mol Syst Biol 2011, 7:495. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Griffith JS: Mathematics of cellular control processes. II. Positive feedback to one gene.
J Theor Biol 1968, 20:209216. PubMed Abstract  Publisher Full Text

Thomas R: Logical analysis of systems comprising feedback loops.
J Theor Biol 1978, 73:631656. PubMed Abstract  Publisher Full Text

Roessler OE: Basic circuits for fluid automata and relaxation systems.

Maruyama T, Li J, Vaque JP, Konkel JE, Wang W, Zhang B, Zhang P, Zamarron BF, Yu D, Wu Y, et al.: Control of the differentiation of regulatory T cells and T(H)17 cells by the DNAbinding inhibitor Id3.
Nat Immunol 2011, 12:8695. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Powell JD, Delgoffe GM: The mammalian target of rapamycin: linking T cell differentiation, function, and metabolism.
Immunity 2010, 33:301311. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Stritesky GL, Muthukrishnan R, Sehra S, Goswami R, Pham D, Travers J, Nguyen ET, Levy DE, Kaplan MH: The transcription factor STAT3 is required for T helper 2 cell development.
Immunity 2011, 34:3949. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Dong C: TH17 cells in development: an updated view of their molecular identity and genetic programming.
Nat Rev Immunol 2008, 8:337348. PubMed Abstract  Publisher Full Text

Feinerman O, Veiga J, Dorfman JR, Germain RN, AltanBonnet G: Variability and robustness in T cell activation from regulated heterogeneity in protein levels.
Science 2008, 321:10811084. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chang HH, Hemberg M, Barahona M, Ingber DE, Huang S: Transcriptomewide noise controls lineage choice in mammalian progenitor cells.
Nature 2008, 453:544547. PubMed Abstract  Publisher Full Text

Wilson HR, Cowan JD: Excitatory and inhibitory interactions in localized populations of model neurons.
Biophys J 1972, 12:124. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Mjolsness E, Sharp DH, Reinitz J: A connectionist model of development.
J Theor Biol 1991, 152:429453. PubMed Abstract  Publisher Full Text

Tyson JJ, Novak B: Functional motifs in biochemical reaction networks.
Annu Rev Phys Chem 2010, 61:219240. PubMed Abstract  Publisher Full Text

Clewley R, Sherwood WE, LaMar MD, Guckenheimer JM:
PyDSTool, a software environment for dynamical systems modeling.
[http://pydstool.sourceforge.net webcite]

Tyson JJ, Novak B: Temporal organization of the cell cycle.
Curr Biol 2008, 18:R759R768. PubMed Abstract  Publisher Full Text  PubMed Central Full Text