Research ArticleImmunology

The Duration of T Cell Stimulation Is a Critical Determinant of Cell Fate and Plasticity

See allHide authors and affiliations

Science Signaling  05 Nov 2013:
Vol. 6, Issue 300, pp. ra97
DOI: 10.1126/scisignal.2004217


Variations in T cell receptor (TCR) signal strength, as indicated by differential activation of downstream signaling pathways, determine the fate of naïve T cells after encounter with antigen. Low-strength signals favor differentiation into regulatory T (Treg) cells containing the transcription factor Foxp3, whereas high-strength signals favor generation of interleukin-2–producing T helper (TH) cells. We constructed a logic circuit model of TCR signaling pathways, a major feature of which is an incoherent feed-forward loop involving both TCR-dependent activation of Foxp3 and its inhibition by mammalian target of rapamycin (mTOR), leading to the transient appearance of Foxp3+ cells under TH cell–generating conditions. Experiments confirmed this behavior and the prediction that the immunosuppressive cytokine TGF-β (transforming growth factor–β) could generate Treg cells even during continued Akt-mTOR signaling. We predicted that sustained mTOR activity could suppress FOXP3 expression upon TGF-β removal, suggesting a possible mechanism for the experimentally observed instability of Foxp3+ cells. Our model predicted, and experiments confirmed, that transient stimulation of cells with high-dose antigen generated TH, Treg, and nonactivated cells in proportions depending on the duration of TCR stimulation. Experimental analysis of cells after antigen removal identified three populations that correlated with these T cell fates. Further analysis of simulations implicated a negative feedback loop involving Foxp3, the phosphatase PTEN, and Akt-mTOR in determining fate. These results suggest that there is a critical time after TCR stimulation during which heterogeneity in the differentiating population of cells leads to increased plasticity of cell fate.


CD4+ T cells can be grouped into two main subtypes: those that mediate the immune response [T helper (TH) cells] and those that suppress the immune response, such as regulatory T (Treg) cells. Treg cells play an important role in suppressing T cell–mediated immunity, limiting the damage caused by the immune response and preventing autoimmune diseases (1); however, Treg cells may play a detrimental role in cancer by acting as an effector of immune suppression by tumors. Patients with cancer often have increased numbers of Treg cells that exhibit greater suppressive function compared to healthy controls, and these features are predictive of poor survival (2). Understanding the factors that control the production of Treg cells has the potential to advance therapies involving either elimination of antigen-specific Treg cells in the context of cancer (3) or enhancement of Treg cell–mediated suppression in the context of autoimmune responses (4).

Treg cells are characterized by the transcription factor forkhead box P3 (Foxp3), a specific pattern of cytokine production, and immunosuppressive function (57). Treg cells suppress other effector T cells through several mechanisms (8), including destruction of effector cells through granzyme B (9); secretion of immunosuppressive cytokines, such as transforming growth factor–β (TGF-β), interleukin-10 (IL-10), and IL-35 (10); metabolic disruption through the production of adenosine (11) or competition for IL-2 (12, 13); and inhibition of dendritic cell (DC) maturation (8, 14). The two main groups of Treg cells that are currently known are natural Treg (nTreg) cells and adaptive (induced) Treg (iTreg) cells. nTreg cells become committed to a regulatory fate while still in the thymus (14), whereas iTreg cells are generated from naïve T cells in the periphery through stimulation by specific factors, including IL-10 (15), TGF-β (16), low-dose antigen (17, 18), and certain DC subsets (1921). Both iTreg cells and TH cells develop in the periphery from common, “naïve,” uncommitted T cell precursors, which differentiate upon encountering cognate peptide-bound major histocompatibility complex (pMHC) on the surface of antigen-presenting cells (APCs).

Previous studies showed that both the concentration of antigen and the duration of stimulation of T cells with antigen markedly influence whether Treg or TH cells develop in response to a particular antigen, such that increased antigen concentrations favor TH cell generation, whereas decreased concentrations favor Treg cell generation (17, 18, 22). Understanding the molecular determinants of this process has major implications for the development of targeted immunomodulation therapies (23). For example, it is important in a cancer vaccination strategy to deliver a large enough dose of antigen to induce the production of tumor-specific TH1 cells. A more delicate balance must be achieved in the context of autoimmunity because many patients have increased numbers of autoreactive T cells, which any therapy must avoid activating. Thus, a strategy designed to induce antigen-specific Treg cells must deliver a dose of antigen that will favor the production of self-antigen–specific Treg cells without activating or inducing autoreactive TH1 or TH17 cells. Few vaccination strategies in cancer or autoimmunity consider antigen concentration, and this may be one reason for the limited success of these efforts to date (24, 25).

The differentiation of Treg cells is influenced by multiple signaling pathways, including those stimulated by engagement of the T cell receptor (TCR), the costimulatory molecule CD28, the IL-2 receptor (IL-2R), and the TGF-β receptor (TGF-βR). Previously, we reported that culturing CD4+ T cells with DCs that present low concentrations of antigen results in both the proliferation of preexisting nTreg cells and the conversion of naïve T cells into iTreg cells (18). Furthermore, the extent of generation of Treg cells inversely correlates with the strength of the TCR signal, as determined by measurement of the phosphorylation of the S6 ribosomal protein downstream of the phosphatidylinositol 3-kinase–Akt–mammalian target of rapamycin (PI3K-Akt-mTOR) signaling pathway, and several studies showed that activation of this pathway inhibits the generation of Treg cells (26, 27). Activation of Akt-mTOR signaling inhibits the development of Treg cells (26), whereas inhibition of the Akt-mTOR pathway enhances Treg cell development (27). Specific components of the Akt-mTOR pathway have also been implicated in the differentiation of TH cell populations, with mTOR complex 1 (mTORC1)–deficient mice failing to generate TH1 and TH17 cells, whereas mTORC2-deficient mice fail to generate TH2 cells (28, 29). In vitro studies demonstrated that the mTORC1 inhibitor rapamycin enables the proliferation of functional Treg cells of murine (30) or human (31, 32) origin.

Mathematical modeling of T cell activation and of signaling pathways in general provides a means to identify important regulatory elements, including feedback loops and feed-forward loops (FFLs), to test hypotheses and to provide testable predictions (33, 34). Modeling TCR activation has so far mainly involved the use of reaction network models that focus on early steps in T cell activation (3538). Because parameterization of reaction network models requires detailed quantitative data that are not currently available for all of the signaling processes involved in T cell differentiation, we have taken an alternative approach based on logical modeling (3941). This approach uses discrete variables to represent network elements, as well as logical rules to describe element interactions. Logical models do not require quantitative parameters, but rather enable the development of complex qualitative networks. Other groups have developed logical models to investigate T cell activation (42, 43) and differentiation (4446), although none of these has explicitly modeled the effects of antigen concentration.

Here, we developed a logical model of T cell differentiation that focuses on both early steps in T cell activation and events leading to the generation or suppression of the Treg cell phenotype. This logical modeling approach enabled us to capture multiple layers of the cellular response to signals, including extracellular stimulation; intracellular signaling that leads to the activation of transcription factors; gene transcription; and autocrine or paracrine responses to secreted factors. Model simulations reproduced known experimental results and made predictions that were confirmed experimentally. Our results suggest that differential activation of signaling pathways downstream of TCR by high and low concentrations of antigen is better modeled as changes in the duration of contact time rather than as changes in the strength of activation of pathways downstream of TCR. Restricting the time of TCR activation in our model resulted in a mixed population of Treg and TH cells, as well as inactive cells, an observation we confirmed experimentally. Analysis of the dynamics of this scenario revealed that the TH cell phenotype stabilized more rapidly than did the Treg cell phenotype, and demonstrated an important role for a feedback loop involving Foxp3, phosphatase and tensin homolog (PTEN), and Akt-mTOR signaling. The logical model presented here demonstrates heterogeneity at the level of the cell population, which recapitulates experimental results and suggests that differential timing of key TCR signaling events plays a major role in the determination of T cell fate.


Naïve T cell differentiation control network

We conducted an extensive literature survey to identify the key signaling events and molecules involved in determining T cell fate upon exposure to antigen. The model that we generated captures critical signaling events, from stimulatory signals by cell-surface receptors to the activation of transcription factors and to the production of proteins critical to specific phenotypes (Fig. 1A). The two primary phenotypes that arise in the model are (i) cells that have the transcription factor Foxp3 but do not produce the cytokine IL-2 and (ii) cells that do not express FOXP3 but have the ability to produce and secrete IL-2. We identify these cells as Treg cells and TH cells, respectively.

Fig. 1 Computational model of T cell differentiation.

(A) Signaling network connections. Pointed arrowheads represent activation, blunt arrowheads represent inhibition, and dashed edges (dashed lines) represent translocation of molecules from the cytosol to the nucleus in response to stimulation of the indicated cell-surface receptors of T cells. (B) Corresponding logic circuit model. Elements are represented by single logic gates (for example, Raf) or combinations of multiple gates (for example, Ras). Element values are computed at each time step on the basis of the values of their inputs with instantaneous updating of gates within elements.

The model includes four receptors: the TCR, the costimulatory receptor CD28, IL-2R, and TGF-βR. The two signals necessary for the activation of naïve T cells are (i) the binding of pMHC complexes on the surface of the APC to the TCR and (ii) the binding of the costimulatory ligands CD80 or CD86 on the APC to CD28 on the T cell (47). The combined signals from the TCR and CD28 activate the transcription factors nuclear factor of activated T cells (NFAT), activator protein 1 (AP-1), and nuclear factor κB (NF-κB). These transcription factors promote expression of IL2 and of the gene encoding the α subunit of IL-2R (IL2RA, also called CD25) (48, 49). The β and γ subunits of IL-2R are constitutively present. Binding of IL-2 to IL-2R results in the activation of another transcription factor, signal transducer and activator of transcription 5 (STAT5) (50, 51), which in turn activates Foxp3, the key transcription factor that drives the differentiation of naïve T cells into Treg cells (52). STAT5 can also modulate the expression of IL2 by inducing production of the transcription factor Blimp-1 (53), but we have not included this interaction in the current version of the model because it complicates the analysis of the simulated dynamics by making IL2 expression transient. We confirmed by additional simulations that including this interaction did not otherwise affect our results.

FOXP3 is expressed after one of the following combinations of transcription factors binds to its promoter: NFAT, AP-1, and STAT5; NFAT and mothers against decapentaplegic homolog 3 (Smad3); or STAT5 alone (54, 55). The Foxp3 protein itself inhibits the expression of IL2, whereas it maintains the expression of IL2RA (7). In parallel with the activation of NFAT, AP-1, NF-κB, and STAT5, the PI3K-Akt-mTOR pathway can also be activated by TCR and CD28 signaling. The activity of the PI3K-Akt-mTOR pathway depends on the activity of PTEN, a phosphatase that is constitutively present in naïve T cells and antagonizes the function of PI3K and is reduced in abundance after TCR activation (56). Reduced PTEN protein abundance or activity results in increased signaling through the PI3K-Akt-mTOR pathway (referred to hereinafter as the mTOR pathway). In addition, PTEN is more abundant in Treg cells than in TH cells (57). On the basis of a report describing the negative impact of mTOR signaling on Treg cell differentiation (58), we assumed that active mTORC1 and mTORC2 can override other signals that stimulate FOXP3 expression. However, because it is well established that TGF-β together with high concentrations of antigen stimulates the generation of Treg cells (16, 59), we also assumed that the combined signals from Smad3, which is activated downstream of TGF-βR, and NFAT stimulate the expression of FOXP3 even when mTOR is active.

Logical model of naïve T cell differentiation

In Boolean network modeling, each node in the graph is represented by a single Boolean variable that takes a value of OFF if the element is inactive and a value of ON if the element is active. For some nodes, however, this two-state representation may be too coarse. For example, because we are studying the effect of varying the TCR signal strength, we need to model at least three possible stimulation values, corresponding to OFF, LOW, and HIGH, to represent the effective antigen concentration. In such cases, depending on the modeling formalism used, node values may be represented by multiple Boolean variables (the approach taken here) or a single discrete variable. We also modeled three values of PI3K to avoid oscillations in the mTOR pathway (see Materials and Methods for additional details).

The interaction map that we generated (Fig. 1A) is not a complete dynamical model because it does not specify how multiple input signals to a node are to be combined. For example, the diagram shows that mTORC2 has two regulators—activation by PI3K and inhibition by ribosomal protein S6 kinase β-1 (S6K1)—but does not specify how the regulations work together. In the logical modeling approach (3941), a dynamical model is constructed with logical functions that determine the next state of each element in the system as a function of its regulators. These functions are constructed from the primitives of Boolean logic—AND, OR, and NOT. For example, mTORC2 activation is described by the following rule: MTORC2′ = PI3K_HIGH or (PI3K_LOW and not S6K1), meaning that mTORC2 is activated either by a HIGH PI3K value that overcomes any inhibition by S6K1 or by a LOW PI3K value in the absence of S6K1.

We defined the full set of elements (table S1) and their corresponding Boolean variables and logic update rules (table S2) with information from our own experiments and from the literature. We tuned the logical functions used in several rules to match the experimental observations (for the references and observations used to determine the regulatory relationships in the model, see table S2). In addition, we set varying thresholds for the elements PTEN, PI3K, mTORC2, and protein kinase Cθ (PKCθ) to match observed differences in the effects of high and low doses of antigen. For example, a high dose of antigen, which produces HIGH TCR activation, is required to inactivate PTEN, which is assumed to be constitutively active in naïve T cells. Because full generation of phosphatidylinositol 3,4,5-trisphosphate (PIP3) and activation of the Akt pathway can only occur if PTEN is inactive, the LOW TCR value that represents a low dose of antigen cannot fully activate either of the mTOR complexes, mTORC1 or mTORC2. As discussed earlier, a HIGH PI3K value is required to overcome the inhibition of mTORC2 activation by S6K1.

We compiled a complete description of the logical model with logic gates that represent the AND and OR rules (Fig. 1B). This representation is similar to the diagrams used in digital circuits and helps to clarify the overall input-output structure of the network. The primary difference between our circuit diagram (Fig. 1B) and those used in digital circuits is that we do not include memory elements for clarity of representation. Memory elements are necessary to keep the state of model elements because all element states are part of the overall system state. In addition to implementing logic functions, logic gates may also have associated delays that we used to model differences in signal propagation between network elements. For example, by including delay gates on the AP-1 pathway, we could fine-tune the timing of transcription factor activation.

We simulated the dynamics of the system with the method of random-order asynchronous Boolean updates (60). Each simulation, or trajectory, represents a sample from the set of possible dynamics that a cell can undergo from a given initial state. Trajectories were started with each variable having a defined value. At each round in the simulation, the next value of each variable was assigned by the logic rules described earlier. Within a round, variables were updated in random order, mimicking the randomness in event timing that may occur because of cell-to-cell variation in the abundance of many proteins and the intrinsic stochasticity of chemical kinetics involving a small number of molecules (61). Time is reported as the number of update rounds, and on the basis of comparison with our experimental results, each round of update corresponds to about 3 hours of real time (see Materials and Methods for additional details). We then summarized the behaviors of the model and how they corresponded to experimental observations (Table 1). The model reproduced most of the experimentally known behaviors of the system and made predictions, several of which we have tested in this work. In addition, an extensive sensitivity analysis (see Materials and Methods for more details) showed that qualitative results did not depend sensitively on the choice of initial conditions or timing parameters (see the “Event timing and delays” section). In the following sections, we discuss these results in more detail.

Table 1 Summary of the observations from the model simulations and data analysis.
View this table:

Architecture of the differentiation control network

The behavior of our model is governed by control loops involving Foxp3. These include positive and negative feedback and coherent and incoherent FFLs, which affect both the transient and steady-state behavior of the model (Fig. 2, A to C). The main FFL begins with signals from the TCR and CD28 propagating either to the IL-2 and CD25 pathway (Fig. 2A, Box 1) or to the mTOR pathway (Fig. 2A, Box 2) and then converging on Foxp3. Because the link from Box 2 to Foxp3 is inhibitory, the signals propagating on parallel branches have opposite signs, leading to designation of this motif as an incoherent FFL (iFFL) (62). A general feature of iFFLs is that they may exhibit transient activation of the output, in this case Foxp3.

Fig. 2 Network architecture and implementation of delay circuits.

(A) Overview of the model control structure. Signaling through TCR and CD28 leads to the production of IL-2 and CD25 (Box 1) and activation of the PI3K-Akt-mTOR signaling pathway (Box 2), which have opposing effects on the production of Foxp3, which itself also exhibits multiple feedbacks. Pointed arrowheads represent activation, blunt arrowheads represent inhibition, and dashed edges represent interactions between boxes. (B) A more detailed view of the control network highlighting several places where delay circuits have been introduced. Pathways leading to the production of IL-2 and CD25 are represented as a sub-box, Box 1.1, and delays in the STAT5-dependent production of Foxp3 and the mTOR-dependent inhibition of Foxp3 activity are represented by d3 and d4, respectively. (C) Details of Box 1.1. The d1 delay replaces the network connections Ras→Raf, Raf→MEK2, MEK2→ERK, and ERK→Fos. The d1′ delay replaces the connections PKCθ→TAK1, TAK1→MKK7, MKK7→JNK, and JNK→Jun. The two delays, d1 and d1′, are always assumed to be equal because activation of AP-1 in the model will always depend on the longer delay on the two pathways. (D) Delay setups (SD1 to SD15) are defined as different combinations of the delay lengths for d1 to d4. The absence of a delay is indicated with an “x,” and the presence of a delay is indicated with a checkmark. For each delay d1(d1′) to d4, we show in the second row the maximum number of rounds that the delay can take, according to the random asynchronous update scheme (see Materials and Methods). (E) Effect of varying delays on the percentage of Foxp3+ cells after stimulation by a high dose of antigen. (F) Effect of varying delays on the percentage of IL-2+ cells after stimulation by a low dose of antigen.

In addition to the iFFL, there exist feedback loops between Foxp3 and each of the FFL branches (Fig. 2, A and B). The mutual inhibition between Foxp3 and mTOR (Fig. 2A) produces a toggle switch (62), which forms the central circuit governing the T cell fate decision. By inhibiting Foxp3 (Fig. 2B), sustained activation of the mTOR pathway leads to a stable TH cell phenotype. On the other hand, activation of Foxp3 activates PTEN, which then inhibits the mTOR pathway, producing a Treg cell phenotype. Positive feedback between Foxp3 and CD25 (Fig. 2A) also stabilizes FOXP3 expression. The inhibition of mTOR by Foxp3 in our model is supported by the experimental observation that TCR-mediated activation of mTOR is attenuated in Foxp3+ Treg cells (18).

As discussed in detail later, the simulation results suggest that the relative timing of several pathways is critical for determining the final state of the toggle switch, and hence the T cell fate decision. These relationships include the time of activation of the mTOR pathway relative to that of the IL-2 and CD25 pathway, as well as the time of activation of Foxp3-PTEN feedback relative to that of signaling from the TCR.

Event timing and delays

The model describes regulatory interactions that can propagate through the network at different speeds. To capture these differences in the logical model, we introduced buffers within key regulatory elements that control the delay between receipt of an input signal and a change in the output state (Fig. 2, B and C, blue boxes). We defined 15 model setups with varying delays, abbreviated SD1 to SD15 (Fig. 2D), and computed time courses after stimulation with high or low doses of antigen (Fig. 2, E and F). Delaying the time at which mTOR inhibits Foxp3 activity (SD3 to SD5, SD8 to SD10, and SD13 to SD15) increases the magnitude of the transient rise in FOXP3 expression (Fig. 2E, top). Fast signal propagation through the AP-1 activation pathway (SD1 to SD5) and the consequent fast expression of IL2 leads to fast activation of STAT5 and expression of FOXP3 before the inhibitory effect of mTOR activation reaches Foxp3, resulting in a large percentage of cells transiently becoming Foxp3+ (Fig. 2E). Conversely, when the delay in AP-1 activation is greater than that of the mTOR pathway (SD11 to SD15, Fig. 2E, bottom, and fig. S1), the percentage of Foxp3+ cells is predicted to decrease.

The model also predicts that cells destined to become Foxp3+ Treg cells transiently become IL-2+ (Fig. 2F). The magnitude of the transient increase in the number of IL-2+ cells is similar for all delay scenarios (Fig. 2F, top), but the timing of the increase can shift to earlier times, depending on the timing of AP-1 activation and its effect on Foxp3 (Fig. 2F, bottom). For the simulations described in the rest of the paper, we used the delay setup defined as SD14 (Fig. 2D) because it best reproduces experimental data for the transient appearance of Foxp3+ cells, as discussed later (see fig. S1 for additional results). Sensitivity analysis (see Materials and Methods for more details) showed that qualitative results did not depend on the choice of delay setup.

Main differentiation scenarios

We validated our model on the basis of its ability to reproduce the main results of several key experiments, namely, those in which cells were treated with a low concentration of antigen alone or a high concentration of antigen alone or in combination with inhibitors or TGF-β. For stimulations with low or high doses of antigen, all trajectories start from a single initial state (table S3) and eventually reach a state from which no further transitions are possible. This state is known as a fixed-point attractor because it represents a point in the state space to which many trajectories converge. Trajectories resulting from stimulation with high and low concentrations of antigen reach distinct fixed points, which are characterized by two sets of element values (Fig. 3A). In addition to the phenotype marker molecules Foxp3 and IL-2, elements that had different values for the two scenarios included PTEN, Akt, and mTORC1, which emphasizes their role in differentiating between the two cell phenotypes. A total of 16 variables had different values and thus distinguished between the high- and low-dose antigen scenarios (Fig. 3A). We selected a subset of 10 key variables (Foxp3, IL-2, PTEN, TCR, Ras, CD25, PI3K, Akt, mTORC1, and mTORC2) that included all of the major pathways that differ between the scenarios, and we removed those variables that were redundant.

Fig. 3 Results for scenarios of high and low doses of antigen.

(A) Values of individual elements in the steady-state attractors that are reached after stimulation with high (H) and low (L) doses of antigen. OFF values are represented by black, whereas ON values are represented by white. (B) Heat map representation of 1000 simulation trajectories for each scenario. (C) Average values of six key elements along the trajectories shown in (B).

The dynamics of the system are shown with trajectory heat maps (Fig. 3B), for which we used a Gray coding scheme (63) to map the 10-element binary vector representing the system state at a given time into a color on the heat map (see table S4 for further details). This approach enables the compact representation of trajectories, in which green and yellow shades represent states with active IL-2, orange and red shades represent states with active Foxp3, and purple and dark blue shades represent states with most elements being inactive. Values for six of the key elements, averaged from the 1000 trajectories computed for each dose, were calculated (Fig. 3C).

Our previous experimental results (18) showed that high concentrations of antigen resulted in the production of mostly TH cells, as demonstrated by the presence of proliferating [cells with low (diluted) amounts of the dye carboxyfluorescein diacetate succinimidyl ester (CFSE)] Foxp3+ cells (Fig. 4A). Our simulations of cell stimulation with a high dose of antigen also showed the predominance of the TH cell phenotype and revealed that a substantial number of cell trajectories passed through a Foxp3+ state (Fig. 4B). We subsequently confirmed in experiments that about 10 to 20% of cells became Foxp3+ 18 hours after stimulation with a high dose of antigen [Fig. 4C; P < 0.05 at 18 hours compared to all other time points; one-way analysis of variance (ANOVA) with Tukey’s posttest]. This information was used to identify the best-fitting delay scenario as SD14 (Fig. 2D, left). Our previous experiments (18) also showed that the abundance of phosphorylated S6 protein (pS6), a surrogate measure of mTOR activity, peaked at 18 hours after stimulation with a high dose of antigen and correlated inversely with the number of Foxp3+ cells remaining after the 7-day culture (Fig. 4C). Our simulation results also recapitulated the negative correlation between pS6 abundance and the percentage of Foxp3+ cells (Fig. 4B).

Fig. 4 Logical model recapitulates the behavior of key elements in the system.

(A) Flow cytometric analysis of a population of CFSE-labeled CD4+ T cells stimulated with a high dose of antigen. Boxes indicate populations of proliferating TH cells and Treg cells. (B) Simulation results for the same conditions as in (A), showing the behavior of Foxp3, pS6, and IL-2 over time. (C) Experimental data showing a transient increase in the percentage of Foxp3+ T cells (squares) after stimulation with a high dose of antigen. The mean fluorescence intensity (MFI) values of pS6 (triangles) are maximal at early time points and are maintained throughout the stimulation (*P < 0.05 at 18 hours compared to all other time points; one-way ANOVA with Tukey’s posttest). (D) Simulation results for Foxp3 for stimulation with a high dose of antigen (HD) alone or with the addition at round 6 of inhibitors of Akt (HD + iAkt), mTOR (HD + imTOR), both Akt and mTOR (HD + iAkt + imTOR), or NFAT (HD + iNFAT). (E) Simulation results for Foxp3, pS6, and IL-2 for stimulation with a high dose of antigen in the presence of TGF-β. (F) Flow cytometric analysis of CD25 and pS6 in gated CD3+ CD4+ T cells stimulated for 18 hours with a high dose of antigen in the absence (left) or presence (right) of TGF-β. Results are representative of three independent experiments, and the numbers in the upper right quadrants represent the means ± SEM of the percentages of CD25+ pS6+ cells from three experiments. (G) Flow cytometric analysis of Foxp3 abundance in a population of gated CD3+ CD4+ T cells stimulated with a low concentration of antigen. (H) Simulation results for the same conditions as in (G) showing the behavior of Foxp3, pS6, and IL-2. (I) Flow cytometric analysis showing the percentages of Foxp3+ Treg cells (squares) and the MFI of pS6 (triangles) in those cells at the indicated time points after stimulation of T cells with a low dose of antigen. Data in (A) and (F) are representative of three independent experiments, whereas data in (C) and (I) are means ± SEM of three independent experiments.

The trajectory heat map for stimulation with a high dose of antigen (Fig. 3B, left, and fig. S2A) revealed three distinct clusters of behavior. Although all trajectories reach the same IL-2+ attractor, about 10% of the trajectories fall into the top cluster, which is characterized by the appearance of red pixels, indicating cells that are transiently Foxp3+. The 80% of trajectories in the middle cluster never become Foxp3+. About 10% of the trajectories fall into the bottom cluster, in which cells are transiently Foxp3+, but then lose both Foxp3 and IL-2 (purple pixels) before becoming IL-2+. The transient appearance of Foxp3+ cells arises from the governing iFFL involving the IL-2 and CD25 pathway and the mTOR pathway, as discussed earlier (Fig. 2A).

A previous study showed that the addition of Akt and mTOR inhibitors to cells 18 hours after they received strong TCR stimulation enhanced the generation of Foxp3+ cells, whereas a similar addition of NFAT inhibitors had no effect (27). As in the experiments, when we modeled inhibition by setting either Akt or mTORC1 to OFF after six rounds of simulation, all of the trajectories produced the Treg cell phenotype (Fig. 4D). Conversely, when the NFAT inhibitor was added after six rounds, there was no increase in the percentage of Foxp3+ cells (Fig. 4D).

TGF-β activates Smad3, which, together with TCR-stimulated NFAT, contributes to the induction of FOXP3 expression by promoting acetylation at the FOXP3 enhancer (64). The activation of Smad3 by TGF-β apparently counteracts the inhibition of Foxp3 by active mTORC1, leading to substantial Treg cell generation, even at high concentrations of antigen (58). Our simulations (Fig. 4E) recapitulated these experimental observations, and they also predicted that Foxp3 abundance increases with faster kinetics in cells that are stimulated with TGF-β together with a high dose of antigen than it does in cells treated with a high dose of antigen alone (Fig. 4C), because the activation of Foxp3 is more direct when TGF-β is added to the stimulation (see Fig. 1). In addition, the model predicts that pS6 remains activated in this scenario because sustained stimulation of the TCR with a high dose of antigen overrides inhibition of mTOR signaling through the Foxp3-PTEN feedback loop (Fig. 2, A and B).

We conducted further experiments to test the prediction that the generation of Treg cells by TGF-β occurs despite sustained mTOR activity, as determined by measurement of pS6 abundance. T cells stimulated with a high concentration of anti-CD3 antibody in the presence or absence of TGF-β exhibited an equivalent abundance of pS6 protein 18 hours after stimulation (Fig. 4F). As expected, TGF-β stimulated a significant increase in the number of Foxp3+ CD4+ cells after 5 days of culture compared to cultures without TGF-β (31.7 ± 10.1% with TGF-β versus 8.8 ± 4.2% without TGF-β; P < 0.01 by one-way ANOVA with Tukey’s multiple comparison posttest). Thus, as predicted by the model, the generation of Treg cells in response to TGF-β does not involve or require inhibition of the mTOR pathway.

Experimental results showed that a low dose of antigen generated an increased percentage of Foxp3+ Treg cells compared to a high dose of antigen (Fig. 4G). Simulation trajectories after stimulation with a low dose of antigen uniformly resulted in the generation of Foxp3+ cells without the generation of pS6 (Fig. 4H). Consistent with this, we observed a reduction in pS6 abundance in cells after 18 hours of stimulation with a low dose of antigen compared to that in cells treated with a high dose of antigen, which was accompanied by an increase in the number of Foxp3+ cells after 7 days (Fig. 4I). On the other hand, whereas 100% of the model trajectories reach the Treg cell steady state (phenotype) at a low dose of antigen, in experiments, only about 35 to 50% of cells were Foxp3+ (Fig. 4G). Thus, the model fails to capture the observed heterogeneity in differentiation outcome. A possible source for this discrepancy is the known variability in the abundance of TCR or other signaling molecules (65) within the naïve T cell population, which is not directly accounted for in our model. Alternatively, variability in the activation timing could also account for the discrepancy.

Varying the duration of TCR stimulation

The failure of the model to reproduce the observed heterogeneity of differentiation outcomes in the situation of stimulation with a low dose of antigen prompted us to explore potential mechanisms beyond changing the activation thresholds between stimulation with low or high dose of antigen. Sauer et al. (27) showed that removal of antigen 18 hours after stimulation resulted in substantial proliferation of Treg cells in culture, which the authors attributed to a reduction in signaling by the mTOR pathway. Other reports have shown that the duration of TCR signaling depends on the concentration of the stimulating antigen (66, 67).

We examined the effect of varying the duration of TCR stimulation by removing the antigen at various times after stimulation, but keeping the same logic rules that were used for stimulation with a high dose of antigen. Under these conditions, simulations reached 11 different fixed-point attractors, and their relative frequencies changed depending on the time of antigen removal (Fig. 5, A and B). These attractors fall into three major cell phenotypes: Treg cell, TH cell, and nonactivated cell (Fig. 5B). If the antigen is removed too early, for example, before time point 2 (T2), only nonactivated cells (IL-2 Foxp3) are generated (Fig. 5A). If the antigen is removed after T10, only TH cells (IL-2+ Foxp3) are observed, whereas antigen removal between T3 and T8 results in increasing numbers of Treg cells (IL-2 Foxp3+), with a peak at T5 and T6 (Fig. 5B and fig. S3).

Fig. 5 Analysis of transient stimulation with a high dose of antigen.

(A) Distribution of attractors when antigen is removed at different times, starting with removal at simulation round 1 (T1), and incrementing the time of removal by one round, up to round 12 (T12). (B) Changes in the sizes of lumped attractors, A1 to A3, A4 to A7, and A8 to A11 (populations characterized as Foxp3 IL-2, IL-2+, and Foxp3+, respectively). (C) Eleven attractors that were observed in the T6 scenario of simulations in which antigen (Ag) was removed are presented in terms of the steady-state values of 10 elements (Foxp3, IL-2, PTEN, TCR, Ras, CD25, PI3K, Akt, mTORC1, and mTORC2) and the number of simulation trajectories ending in each given attractor for 1000 simulation runs. Colors correspond to those in (A). (D) Flow cytometry plot identifying the three major populations of cells produced after activation. Purified CD4+ T cells were stimulated on plates coated with anti-CD3 antibody (1 μg/ml) in the presence of soluble anti-CD28 antibody. After 18 hours, the cells were removed from the coated plates, placed in new wells, and cultured for a total of 96 hours. Cells were analyzed by flow cytometry, and the results show the detection of CD25 and Foxp3 in gated CD3+ CD4+ cells. CD25 Foxp3 cells (purple box) are nonactivated T cells; CD25+ Foxp3 cells (green box) are TH cells; and CD25+ Foxp3+ cells (orange box) are Treg cells. (E) Purified CD4+ T cells were stimulated as described in (D), and the cells were removed from stimulation at the indicated times. All cells were cultured for a total of 96 hours and then were analyzed by flow cytometry to identify nonactivated cells (No Act, purple circles), TH cells (green squares), and Treg cells (orange triangles), as defined in (D). Data in (D) are representative of three independent experiments. Results in (E) are means ± SEM of four to eight replicates from three independent experiments.

We selected the T6 time point for further analysis because it generates all 11 attractors (A1 to A11), which can be distinguished by variations in the state of key element values (Fig. 5C). Attractors A1 to A3 correspond to nonactivated cells (Fig. 5C), whereas attractors A4 to A7 and A8 to A11 represent the TH cell and Treg cell phenotypes, respectively. Of the TH cell attractors, A5 is the most common outcome and is nearly identical to the attractor for a high dose of antigen (Fig. 5C). The Treg attractors A9 to A11 occur with equal frequency, and A11 is most closely related to the attractor for a low dose of antigen (Fig. 5C). The attractors A9 and A10 differ from A11 and the attractor for a low dose of antigen in the elements Ras (ON in A10 and A11, OFF in A9) and PI3K and mTORC2 (ON in A11 and OFF in A9 and A10). These results suggest that the state values of Ras, PI3K, and mTORC2 are not critical to the outcome, whereas the state values of Akt and mTORC1 are more important for the Treg cell phenotype.

The presence of 11 different cell fates after the removal of antigen in our model prompted us to determine whether these could be identified experimentally. Thus, we activated purified CD4+ T cells with a high concentration of immobilized anti-CD3 antibody and removed the cells from the stimulation at various time points after activation. When cells were removed from antigen after 18 hours and then were analyzed after an additional 78 hours of culture without stimulation, we observed three distinct populations of cells (Fig. 5D). Consistent with the simulations, we observed CD25 Foxp3, CD25+ Foxp3, and CD25+ Foxp3+ cells, representing nonactivated cells, TH cells, and Treg cells, respectively. Varying the duration of antigen stimulation from 3 to 96 hours resulted in a distribution of phenotypes (Fig. 5E) in agreement with the modeling results (Fig. 5B).

Dynamics of the cell fate decision

The observation that the removal of antigen at T6 results in 11 different possible outcomes gave us the opportunity to analyze in detail the factors that affect cell fate selection. We generated a heat map of 1000 simulation trajectories for this scenario (Fig. 6A), which are clustered by steady-state attractors sorted from A1 (bottom rows) to A11 (top rows). We also generated additional heat maps showing system trajectories in scenario T6 for different delay setups (fig. S4B). Differentiation into the TH cell phenotype appears to occur faster than does differentiation into the Treg phenotype; thus, we calculated the mean times required to reach steady state. We found that, on average, it took eight simulation rounds longer to reach the Treg cell (Foxp3+) phenotype than it did to reach the TH cell (IL-2+) phenotype (Fig. 6B). This held true for all attractors that led to IL-2+, rather than Foxp3+, cell fate (Fig. 6B).

Fig. 6 Analysis of the timing effects in the scenario of antigen removal.

(A) Trajectories of 1000 model simulations for the scenario T6 (antigen removal at round 6) showing different transient behaviors and different steady states. (B) Average numbers of simulation rounds required to reach steady state for the dominant attractors for scenario T6 and for high and low doses of antigen (error bars show ± SEM). (C) Left: Average values for six key elements leading to attractor A5 (top) and attractor A11 (bottom) computed from 374 and 175 trajectories, respectively. Right: Diagrams summarizing the order of events along trajectories leading to attractor A5 (top) and attractor A11 (bottom). Corresponding plots and diagrams for the remaining attractors are shown in fig. S2. (D) Values of six elements from (C) averaged across simulations leading to each of the attractors (374 for A5, 1000 for high dose, and 175 for A11) for scenario T6, in rounds 3, 4, and 5 (before antigen removal), round 6 (when antigen is removed), round 7 (transient round), and round 18 (when all trajectories reach steady state).

The differences in the speed of reaching the Foxp3+ and IL-2+ attractors arise from differences in the order in which feedback loops and FFLs are activated (Fig. 6C, right), which in turn is a result of the stochasticity inherent in our simulations. IL-2+ phenotypes result from situations in which the effects of antigen removal propagate slowly enough to permanently shut off FOXP3 expression (Fig. 6C, top). Foxp3+ phenotypes have two main requirements: (i) early (transient) activation of Foxp3 that occurs through the CD25-STAT5 pathway and (ii) early activation of Foxp3 that leads to activation of PTEN and inhibition of mTOR. Foxp3 IL-2 phenotypes (figs. S5 and S6) arise from either (i) slow propagation of the TCR ON signal or (ii) rapid inhibition of activation after antigen removal.

Our analysis shows that the critical factor distinguishing between trajectories leading to A5 (the dominant IL-2+ attractor) and A11 (the dominant Foxp3+ attractor) is transient induction of FOXP3 expression in the rounds just after antigen removal (Fig. 6, C and D). Although the fraction of Foxp3+ trajectories is less than 100% for any specific round, analysis of the individual trajectories finds that all trajectories leading to A11 are transiently Foxp3+ at some point after antigen removal. Foxp3 begins to activate PTEN at round 8, and this activation is complete by round 11. After this, PTEN remains active because of the absence of further TCR stimulation, and the activity of mTORC1 decreases rapidly, which enables expression of FOXP3 to be resumed, which stabilizes the Treg cell phenotype (Fig. 6C).

We wondered whether we could predict the subsequent behaviors of the A5 and A11 attractors from the state of the cell at the time of antigen removal (Fig. 6D). The increase in mTORC1 activation appeared to occur more slowly in the A11 trajectories than in the A5 trajectories, such that by round 6, the average mTORC1 activity was lower in A11 than in A5 or in trajectories for a high dose of antigen (Fig. 6D). We measured the extent of mTORC1 activity experimentally by assessing the percentages of cells containing pS6, which is downstream of mTORC1 activation. After 6 hours of activation of CD4+ T cells with a high concentration of anti-CD3 antibody, we observed a substantial number of pS6+ cells, but no CD25+ cells (Fig. 7A, left). After 18 hours of activation, we detected CD25+ cells, some of which were also pS6+, whereas others were pS6 (Fig. 7A, right). We determined the relative percentages of pS6+, CD25+, and pS6+ CD25+ cells over time and observed an early peak of cells that were pS6+ in the absence of CD25 (Fig. 7B). There was a peak in the percentage of pS6 CD25+ cells between 12 and 18 hours after activation, and the number of pS6+ CD25+ cells increased progressively after 12 hours (Fig. 7B). We analyzed the same populations of cells in the model and observed very similar kinetics for these three populations in the simulations (Fig. 7C).

Fig. 7 State of cells at the time of antigen removal.

(A) Flow cytometric analysis of the activation state of T cells stimulated for 6 and 18 hours. Purified CD4+ T cells were stimulated on plates coated with anti-CD3 antibody (1 μg/ml) in the presence of soluble anti-CD28 antibody. After 6 or 18 hours, the cells were analyzed by flow cytometry, and cells gated on CD3 and CD4 were analyzed for the presence of pS6 and CD25. (B) Time course of the activation state of T cells stimulated with a high concentration (1 to 2 μg/ml) of anti-CD3 antibody in the presence of constant anti-CD28. Cells were analyzed at the indicated times for the presence of pS6 and CD25. Three populations of cells were followed over time: pS6+ CD25 (black circles), pS6+ CD25+ (brown squares), and pS6 CD25+ (orange triangles). Data are means ± SEM of four to eight replicate samples from three independent experiments. (C) Distribution of pS6+ CD25, pS6+ CD25+, and pS6 CD25 simulation trajectories for stimulation with a high dose of antigen at different rounds before the antigen was removed.


Here, we developed a logical network model of TCR signaling (Fig. 1) to investigate the role of antigen concentration in determining T cell fate after the activation of critical signaling pathways downstream of TCR, CD28, IL-2R, and TGF-βR. The model reproduced important experimental results (Table 1), including the generation of Foxp3 TH cells in response to stimulation with a high dose of antigen, as well as the generation of Foxp3+ Treg cells in response to a low dose of antigen (Fig. 3). In addition, our model captured the experimentally observed inverse correlation between the generation of Treg cells and the activation of the mTOR pathway, as determined by measurement of S6 phosphorylation (Fig. 4, B, C, H, and I). The model also made several predictions that we verified experimentally (Table 1), including the transient appearance of Foxp3+ cells within a population of T cells activated with a high dose of antigen (Fig. 4, B and C) and the sustained activation of mTOR in Treg cells in response to a combination of TGF-β and a high dose of antigen (Fig. 4, E and F). The model also predicted, and experiments confirmed, that transient stimulation of CD4+ cells with a high dose of antigen would result in three stable T cell fates (TH cells, Treg cells, and nonactivated cells), and that the relative proportions of these cells depended on the duration of stimulation (Fig. 5). Experimental analysis of the cell population at the time of antigen removal identified three distinct populations on the basis of the cell-surface abundance of CD25 and the extent of mTOR activation, which correlated with these T cell fates (Fig. 7). Overall, our model provides insights into how TCR signal strength influences T cell differentiation and demonstrates an important role for the differential timing of key signaling elements in determining cell fate.

We observed in simulations that a proportion of T cells exposed to a high dose of antigen transiently become Foxp3+ (Fig. 4B). An early increase in the number of Foxp3+ cells after stimulation with a high dose of antigen was shown previously for human (68) and mouse (69) T cells. This finding can be explained by an iFFL starting with TCR activation and passing through the activation of the IL-2–CD25–STAT5 and PI3K-Akt-mTOR pathways (Fig. 2A). Upon stimulation of CD4+ T cells with a high dose of antigen, FOXP3 can be transiently induced by STAT5, but it is eventually blocked by activation of the mTOR pathway. We tested this prediction from the model in our experimental system with mouse T cells, and we observed that a proportion of T cells became transiently Foxp3+ upon stimulation with a high dose of antigen (Fig. 4C). A more detailed analysis of these activation-induced Foxp3+ T cells revealed that these cells do not have suppressor function and that they exhibit increased methylation of the Foxp3 locus (69); however, our simulations also revealed that the transient increase in the number of Foxp3+ cells plays a critical role in the differentiation of Treg cells in situations of low or unsustained TCR signaling.

Our model also suggests a possible explanation for the apparent instability of Treg cells generated by TGF-β. It is well established that TGF-β stimulates the production of Treg cells when used in combination with a high dose of antigen (16, 20). Treg cells generated in this fashion are Foxp3+ and suppressive, but data from another study suggest that this phenotype is not stable after removal of the TGF-β, and that it is associated with incomplete demethylation of the Foxp3 locus (70). Simulations of the model in which TGF-β is present together with a high dose of antigen show the rapid generation of Foxp3+ cells, which is accompanied by sustained activation of the mTOR pathway (Fig. 4E). This simulation result is in contrast to our findings from experiments and model simulations with a low dose of antigen in the absence of TGF-β, which showed that the generation of Foxp3+ cells was associated with reduced mTORC1 activation (Fig. 4H). We tested this prediction experimentally and observed that the proportion of cells with active pS6, a marker of mTORC1 activity, was maintained in the presence of TGF-β, despite the generation of substantial numbers of Treg cells in these cultures (Fig. 4F). These results suggest that the addition of TGF-β overrides the negative effects of mTOR signaling on FOXP3 expression, and provide a potential mechanism by which depletion of TGF-β may lead to the suppression of FOXP3 expression by activated mTORC1. The continued activation of the mTOR pathway in this setting may thus provide an explanation for the instability of Treg cells generated by TGF-β under conditions of enhanced TCR signaling.

The logical modeling approach enabled us to model complex behaviors with a simple approach, with two values for each variable. Several logical models of T cell differentiation and activation were developed previously (4346, 71). Two of these models (43, 71) considered TCR signaling in detail, but did not examine the ultimate fate of the activated T cell. Other models describing TH cell differentiation (4446) have focused on the interplay between cytokines that stimulate the differentiation of specific TH cell subsets, but they did not treat molecular pathways in detail, nor did they analyze the effect of TCR signal strength. Thus, our model bridges a gap in the existing literature by connecting TCR signaling cascades to the fate of the differentiating T cell. A limitation of our model is that it considers only two major outcomes, Treg cells and TH cells, without tracking the further differentiation of TH cells into the TH1, TH2, or TH17 subset. It is known that the concentration of antigen has a profound influence on the differentiation of TH1 and TH2 cells (66), and we plan to further develop the existing model to include these cell fates in the future.

Although all elements of the model were initially given only two values (ON or OFF), we observed that this could give rise to oscillatory behavior in the scenario of a low dose of antigen in which a steady state is never reached. We overcame this limitation by enabling two of the elements (PI3K and PIP3) to have three values (OFF, LOW, and HIGH) and by setting different activation thresholds for the low and high doses of antigen. These thresholds were centered on the ability of the TCR signal to inhibit the critical negative regulator PTEN and to activate PKCθ as well as the mTOR pathway. These initial attempts to model the generation of Treg cells after stimulation with a low dose of antigen resulted in a model that was too rigid, producing 100% Foxp3+ cells (Fig. 4H), instead of the mixed population of 35 to 50% Treg cells that we observed experimentally (Fig. 4G).

We subsequently used an alternative interpretation of stimulation with a low dose of antigen in which the duration of TCR activation was varied while the same activation threshold was maintained. Sauer et al. (27) showed that the removal of antigen after a short period of stimulation resulted in the substantial proliferation of Treg cells in culture, which was attributed to a reduction in mTOR signaling. By recapitulating this scenario with activation thresholds for a high dose of antigen, but varying the length of time that the TCR signal was on, we found an optimal duration of TCR signaling that generated both TH cell and Treg cell phenotypes in about the same proportions that we observed experimentally (Fig. 5, A and B). We confirmed this prediction experimentally by showing that removing T cells from the activating stimulus at various times also generated the same pattern of T cell differentiation (Fig. 5, D and E). This result suggests that differential activation of signaling pathways downstream of TCR by high or low doses of antigen can arise from changes in the duration of contact time between T cells and APCs, which is consistent with several other studies (22, 66). These reports showed that the generation of Treg cells in vivo by a low dose of antigen requires a high-affinity pMHC-TCR interaction, whereas low-affinity pMHC-TCR interactions are ineffective (22, 66). This difference is related to the induction of IL-2 responsiveness, such that, in contrast to low concentrations of high-affinity peptides, high concentrations of low-affinity peptides failed to stimulate substantial increases in the cell-surface abundance of CD25, and thus cells activated in this way exhibited reduced STAT5 signaling in response to IL-2 (66). Imaging studies also revealed that only a small percentage of T cells that sense low-affinity pMHC have long-lasting, stable contacts with DCs, compared to those of T cells that sense low concentrations of high-affinity pMHC (66). Another study showed that a single antigen-specific T cell in vivo can differentiate into different TH cell subsets, and that this property is dependent on the concentration of the antigen delivered to the cell (72).

Analysis of individual simulation trajectories enabled us to identify critical elements that contribute to the determination of T cell fate (Figs. 6 and 7). Two main findings emerged from this analysis. The first was the importance of the feedback loop (Fig. 6C) involving IL-2R, PTEN, Foxp3, and Akt-mTOR signaling. Trajectories leading to the TH cell phenotype have increased mTORC1 activation before the antigen is removed, thus enabling inhibition of Foxp3 (Fig. 6D, A5 trajectories). Cells that become Treg cells have decreased mTORC1 activity, which enables FOXP3 expression to be induced and leads to the increased activity of PTEN. At later times, PTEN activation leads to inhibition of mTORC1 signaling, enabling the sustained generation of Foxp3+ cells. Experimentally, we showed that mTORC1 activation occurred very early, before the appearance of cell-surface CD25, and that at later time points, cultures became dominated by CD25+ cells with increased mTORC1 activation. At intermediate time points, we identified an additional cell type that was CD25+ without mTORC1 activation (Fig. 7B). We hypothesize that the heterogeneity in T cell activation status observed at these intermediate time points contributes to the plasticity in cell fate that we observed when activation was interrupted at these times.

The second major finding from the trajectory analysis is that the differentiation of cells into Treg cells takes more time than does differentiation into TH cells (Figs. 6B and 7). This suggests that there is fundamentally more instability in the pathway toward Treg cell differentiation, and that TH cells switch more or less directly from the naïve state to the final TH phenotype, whereas cells destined to be Treg cells pass through an additional intermediate state in which they produce IL-2 before PTEN is activated, which enables sustained activation of Foxp3.

PTEN, which emerged as a critical regulator in this analysis, plays a critical role in modulating the activity of the mTOR pathway, but the regulation of the abundance and activity of PTEN is not well understood. Two potential mechanisms for the regulation of PTEN are included in the current version of the model: A high concentration of antigen decreases the activity of PTEN, whereas Foxp3 activity increases the activity of PTEN. These mechanisms are consistent with experimental observations that Treg cells have more PTEN than do other T cells, and that Treg cells signal poorly through the Akt-mTOR pathway (73). Although we are not aware of direct evidence that Foxp3 causes an increase in abundance of PTEN, previous work showed that PTEN is responsible for the reduced activation of mTOR in Treg cells (73), and that the activation of Treg cells causes a substantial increase in PTEN abundance. The combination of transient FOXP3 expression with the ability of Foxp3 to increase PTEN activity was responsible for the generation of Treg cells in scenarios in which antigen was removed, suggesting one possible mechanism for the stabilization of the Foxp3+ phenotype. Further experimental work is necessary to define the mechanisms by which Foxp3 leads to increases in PTEN abundance and activity.

In conclusion, our findings suggest that logical modeling can provide insights into the role of TCR signal strength in determining the outcome of T cell differentiation. The model made several predictions that we confirmed experimentally, and further refinements to the model will increase its power and usefulness to immunologists. In the future, this model can be refined to aid in the design and development of targeted immunotherapies for diseases, such as autoimmunity, cancer, and transplant rejection, in which the TH cell subsets that are generated play a critical role in disease outcome.

Materials and Methods

Discrete logical model development

To study the control of T cell differentiation, we created a logic circuit model following an iterative procedure (74). First, we defined key elements in the system to include in the model. The next step was to define element regulatory sets, as well as the number of values that were necessary to describe each element. As we discussed earlier, we aimed to restrict the number of values to the minimum that still enabled us to reproduce experimental results. Therefore, we initially chose all discrete variables that represent individual elements to have two values (OFF and ON), except for the TCR variable, which needs to have three values (OFF, LOW, and HIGH) to implement scenarios of low and high doses of antigen. Because the negative loop between mTORC1 and mTORC2 includes elements with only two values, it always oscillates. To avoid oscillations that do not occur experimentally, we decided to model PI3K (and consequently, PIP3), which is upstream of the mTORC1-mTORC2 loop (and, in our model, is a direct regulator of mTORC2 and an indirect regulator of mTORC1), with three values, similar to the TCR. Another solution would be to model all of the elements in the loop with three values, as we plan to do in the future. For those elements that are implemented with multivalued variables, translation of a discrete variable into a vector of Boolean variables is necessary. In our case, this is simple, because it includes translating from three-valued to two-valued variables, that is, defining the LOW and HIGH variables for TCR, PI3K, and PIP3. When developing logic rules, we decided how to include LOW versus HIGH variables in the rules according to our experimental results and a review of the literature.

Model simulation method

Having built the model, we ran simulations to determine whether it could accurately capture the experimental results that we obtained previously. Simulations were performed with the BooleanNet tool (60). Logical rules were updated in each simulation step using a random asynchronous approach, where in each step, only one variable was chosen and the rule associated with that variable was evaluated to obtain a new value for the variable. The choice of rule to be evaluated next can follow different algorithms. The procedure that is applied in the BooleanNet tool is a random choice of the order of rules in one update round. That is, in each update round, all rules are accessed and evaluated once, but the order in which this is done is random. We simulated and analyzed several scenarios, and for each scenario, 1000 independent simulations were conducted. Each run consisted of up to 30 update rounds. The stable presence of Foxp3 was used as a marker for the Treg cell phenotype, and the stable presence of IL-2 was used as a marker for the TH cell phenotype. To run simulations, it was first necessary to decide the initial conditions, that is, to determine the initial variable values when the simulation starts. Because we were interested in studying the differentiation of naïve T cells and understanding the rationale behind the differentiation, we started with the variable values that best represented the naïve T cell state. Most of the variables were assumed to initially have a logical “0” value, which means that they were inactive or absent in the beginning. The only elements and their corresponding variables that we assumed to be active or present initially were TCR (either LOW or HIGH variable, depending on the scenario we were analyzing), CD28 (stimulation without costimulation is suggested to lead to anergy), CD122 and CD132 (assumed to be constitutively active), PTEN (suggested to have a nonzero value in naïve T cells), and TSC1-TSC2 (assumed to be ON, such that it inhibits Rheb until Akt is activated by TCR signaling).

Analysis of simulation trajectories

To emphasize the stochastic behavior of the system leading to different trajectories for the same initial conditions and same stimulation scenario, we visualized the computed system trajectories with heat maps (Figs. 3B and 6A). Each color in a heat map represents a different value of a vector of 10 Boolean variables implementing key elements of the model (Foxp3, IL-2, PTEN, TCR, Ras, CD25, PI3K, Akt, mTORC1, and mTORC2). We selected these elements on the basis that they can be used to characterize the system at any time point with minimal loss of information. A heat map shows 1000 simulation trajectories of the same scenario (same stimulation and initial state of the 10-element vector), with each row representing a single trajectory (changes in the 10-element vector from initial state to simulation round 30). To group similar trajectories in the map, we assumed that values of the 10-element binary vector encoded decimal values with Gray code, a binary number system in which two successive values differ in only one bit (see table S4). In contrast to binary code, Gray code is a nonweighted code. For example, the four-digit binary numbers 0111 and 1000 in binary code encode the decimal numbers 7 and 8, whereas the same decimal numbers are encoded by 0100 and 1100 in Gray code. Computation of decimal values according to the Gray code enabled us to group simulation trajectories that were closest in terms of individual element behaviors because two consecutive decimal numbers will always represent vectors that differ in a single element. The grouping of trajectories was performed by sorting trajectories according to decimal vector values starting from the end of the trajectory (the last simulated state or steady state) toward the beginning (the initial state). This sorting approach enabled the grouping of all trajectories leading to same steady state together (Fig. 6A), and then, within these groups, to cluster closer those trajectories that have similar behavior. To emphasize differences across studied scenarios, we studied individual element behavior by averaging element values in each round across 1000 trajectories, from initial state to round 30. We then compared the average trajectories of different elements within the same scenario, as well as the trajectories of the same element across different scenarios.

Sensitivity analysis

We conducted an extensive analysis of the sensitivity of our model, in the scenarios described earlier, to changes in internal model delays, to turning individual model elements OFF or ON, and to different initial element values. Changes in model delays can be interpreted as variations of model parameters. As expected, such changes lead only to quantitative differences in the average values of individual model elements, but not to qualitative differences in overall system behavior and observed phenotypes (Fig. 2, E and F, and figs. S1, S2, and S4). On the other hand, similar to our study of inhibitor addition at simulation round 6, switching elements of the model permanently OFF affects logical rules and can lead to qualitative differences in observed phenotypes (fig. S7). We expanded our study of the scenario of inhibitor addition in the context of a high dose of antigen and showed that timing did not play an important role in the case of turning model elements OFF (fig. S7E). Our assumption of random initial values for model elements did not affect the steady-state values in the scenarios of a high dose of antigen alone or in the presence of TGF-β (fig. S8, first and second columns) but resulted in a mixed population for stimulation with a low dose of antigen (fig. S8, third column) and led to a small percentage of activated cells in the absence of stimulation (fig. S8, fourth column).

Purification and activation of DCs and CD4+ T cells

DCs were generated as previously described (18) by growing bone marrow precursors for 5 days in RPMI medium containing granulocyte-macrophage colony-stimulating factor (5 ng/ml) and IL-4 (5 ng/ml). DCs were purified with a 13.5% Histodenz gradient (Sigma), as previously described, to enrich for the mature population. CD4+ T cells were purified from the dissociated spleens and lymph nodes of BDC2.5 mice by negative selection with a CD4+ T cell isolation kit (Miltenyi Biotec). DCs and T cells were cocultured in 96-well U-bottomed plates (BD Falcon). DCs were seeded at 2 × 104 cells per well, and peptides were added at least 2 hours before the addition of CFSE-labeled CD4+ T cells (1 × 105 per well). Cells stimulated with AV10 peptide at a concentration of 0.01 ng/ml were considered to be exposed to a low dose of antigen, whereas cells stimulated with AV10 peptide at a concentration of 100 ng/ml were considered to be exposed to a high dose of antigen. Cultures were maintained in complete Dulbecco’s modified Eagle’s medium containing 10% fetal bovine serum (FBS), penicillin and streptomycin, l-glutamine, sodium pyruvate, nonessential amino acids, Hepes buffer, and β-mercaptoethanol. Cells were fed on days 3, 4, 5, and 6 by replacing half of the culture medium with fresh medium. Cells were harvested at the times indicated in the figure legends for flow cytometric analysis. In experiments with TGF-β or in which antigen was removed, purified CD4+ T cells were stimulated on tissue culture wells precoated with anti-CD3 antibody (1 μg/ml) and with soluble anti-CD28 antibody in the presence or absence of TGF-β (5 ng/ml). In experiments in which antigen was removed, cells were removed from the coated plate at various time points and were placed in fresh uncoated plates for the remainder of the culture period. Cells were analyzed by flow cytometry for the presence of pS6 at various time points and for Foxp3 4 to 5 days after stimulation.

Cell staining and flow cytometry

Cell-surface staining was performed on cells resuspended in phosphate-buffered saline containing 2% FBS with the following antibodies: phycoerythrin-conjugated anti-CD11c antibody, allophycocyanin-conjugated anti-CD25 antibody (eBioscience), and peridinin chlorophyll protein–conjugated anti-CD4 antibody (BD Biosciences). Intracellular Foxp3 protein was detected with a Pacific Blue–conjugated anti-Foxp3 antibody and staining buffers (eBioscience) according to the manufacturer’s instructions. Staining with antibodies against pS6 (Cell Signaling Technology) was performed in conjunction with the Foxp3 antibody and staining buffers. Stained cells were analyzed on an LSR II flow cytometer (BD Biosciences) with FACSDiva software (BD Biosciences).

Statistical analysis

Statistics were performed with GraphPad Prism software (GraphPad Software) using one-way ANOVA with Tukey-Kramer multiple comparison posttest. Results were considered significant if P < 0.05.

Supplementary Materials

Fig. S1. Behaviors of Foxp3 and IL-2 for different delay setups in scenarios of high and low doses of antigen.

Fig. S2. System trajectories for different delay setups in scenarios of high and low doses of antigen.

Fig. S3. Frequency with which attractors A1 to A11 occur in the antigen removal scenarios T1 to T12.

Fig. S4. Analysis of the effects of different delay setups for the scenario of stimulation with a high dose of antigen followed by antigen removal.

Fig. S5. Signal propagation and element trajectories of different attractors in scenario T6.

Fig. S6. Transient behavior and phenotypes in the presence of inhibitors and a high dose of antigen.

Fig. S7. Analysis of the effects of switching model elements OFF.

Fig. S8. Analysis of the effects of randomizing initial values of model elements.

Table S1. Model elements.

Table S2. Model rules.

Table S3. Initial element values for several scenarios.

Table S4. Computation of numbers for heat maps.

References for tables S1 and S2.

References and Notes

Acknowledgments: We thank J. Sekar for helpful conversations in the construction of the model. Funding: This work was supported by NIH grants P01CA73743 (P.A.M.), UL1-RR024153 (N.M.-Z. and J.R.F.), and GM080398 (L.P.K.); American Diabetes Association grant 1-06-RA-94 (P.A.M.); and National Science Foundation Expeditions in Computing grant 0926181 (N.M.-Z. and J.R.F.). M.S.T. was supported by NIH training grant 5T32 CA82084-10. Author contributions: N.M.-Z., P.A.M., and J.R.F. conceived and designed the experiments; M.S.T. and P.A.M. performed the experiments; N.M.-Z. and J.R.F. developed and tested the logical model; N.M.-Z., P.A.M., and J.R.F. analyzed the data; L.P.K. contributed reagents, materials, and analysis tools; and N.M.-Z., P.A.M., and J.R.F. wrote the paper. Competing interests: The authors declare that they have no competing interests.
View Abstract

Stay Connected to Science Signaling

Navigate This Article