A Mechanistic Computational Model Reveals That Plasticity of CD4+ T Cell Differentiation Is a Function of Cytokine Composition and Dosage

CD4⁺ T cells provide cell-mediated immunity in response to various antigens. During an immune response, naïve CD4⁺ T cells differentiate into specialized effector T helper (Th1, Th2, and Th17) cells and induced regulatory (iTreg) cells based on a cytokine milieu. In recent studies, complex phenotypes resembling more than one classical T cell lineage have been experimentally observed. Herein, we sought to characterize the capacity of T cell differentiation in response to the complex extracellular environment. We constructed a comprehensive mechanistic (logical) computational model of the signal transduction that regulates T cell differentiation. The model’s dynamics were characterized and analyzed under 511 different environmental conditions. Under these conditions, the model predicted the classical as well as the novel complex (mixed) T cell phenotypes that can co-express transcription factors (TFs) related to multiple differentiated T cell lineages. Analyses of the model suggest that the lineage decision is regulated by both compositions and dosage of signals that constitute the extracellular environment. In this regard, we first characterized the specific patterns of extracellular environments that result in novel T cell phenotypes. Next, we predicted the inputs that can regulate the transition between the canonical and complex T cell phenotypes in a dose-dependent manner. Finally, we predicted the optimal levels of inputs that can simultaneously maximize the activity of multiple lineage-specifying TFs and that can drive a phenotype toward one of the co-expressed TFs. In conclusion, our study provides new insights into the plasticity of CD4⁺ T cell differentiation, and also acts as a tool to design testable hypotheses for the generation of complex T cell phenotypes by various input combinations and dosages.

Introduction

The diversity and number of immunity-related diseases require a high level of heterogeneity in the immune system to maintain the overall well-being of a human. Early studies of immune responses led to a discovery that the CD4⁺ T cells (referred to as T cells), which are critical players in immunity, can be classified into two subtypes - T helper 1 (Th1) and T helper 2 (Th2) cells (Mosmann et al., 1986). Each type of effector T cell produces a specific set of cytokines that define the function of the cell and the way it further governs the immune response. Specifically, the Th1 cells are responsible for several autoimmune diseases, whereas the Th2 cells are the mediators in cases of allergy and asthma (Reiner, 2007; Zhu and Paul, 2008). More recently, a number of additional T cell subtypes, including the inducible regulatory T cells (iTregs) (Groux et al., 1997; Chen et al., 2003; Schmitt and Williams, 2013), T helper 17 (Th17) (Romagnani, 2000; Harrington et al., 2005; Mangan et al., 2006), T helper 9 (Th9) (Dardalhon et al., 2008; Veldhoen et al., 2008; Soroosh and Doherty, 2009), and follicular T helper cells (Tfh) (Breitfeld et al., 2000; Schaerli et al., 2000) have been discovered, and their functions have been extensively studied. For example, the Th17 cells have been found to be responsible for assisting the immune response against extracellular bacteria and fungi, whereas the main role of the iTregs is to maintain the balance and regulate immune responses by the T helper cell subtypes (Zhu and Paul, 2008). The Th9 cells have been found to be involved in pathogen immunity and inflammatory diseases (Kaplan, 2013). Finally, the Tfh cells assist in T cell-dependent B cell response (Breitfeld et al., 2000; Schaerli et al., 2000; Ma et al., 2012).

In addition, recent studies suggest that some T helper cells are capable of switching and exhibiting phenotypes of one of the alternative effector T cells, depending on the combination of input signals that the cell receives. For example, the iTregs and Th17 can switch from one phenotype to the other in response to the pleiotropic cytokine interleukin-6 (IL-6) (Xu et al., 2007; Lee et al., 2009a; Rowell and Wilson, 2009; Kimura and Kishimoto, 2010). The fully differentiated Th17 cells have been observed to produce Th1-cell-specific cytokines (Shi et al., 2008; Lee et al., 2009b; Nindl et al., 2012; Harbour et al., 2015). The Th2 cells have been reported to further develop into Th9 cells (Veldhoen et al., 2008). More complexity in T cell differentiation was observed in the form of co-expression of mutually exclusive lineage-specifying transcription factors (TFs) (Peine et al., 2013; Bock et al., 2017). This co-expression can lead to the development of stable or intermediate subtypes that share characteristics of more than one type of T cell (Tartar et al., 2010). Examples of such mixed (complex) phenotypes include Th1–Th2 (Peine et al., 2013; Bock et al., 2017) and Th1–Th17 (Kullberg et al., 2006; Morrison et al., 2013).

The differentiation process is governed by the regulation of multiple, mutually cross-linked signaling pathways, which form complex networks (Zhu et al., 2010). The stimulation of the naive CD4⁺ T cells by various cytokines triggers a cascade of signaling events, such as the activation of the JAK/STAT pathways that lead to the activation of T cell lineage-specifying TFs (Murphy and Reiner, 2002; Kaiko et al., 2008). For example, the commitment to Th1 lineage is initiated through signaling by interferon gamma (IFN-γ) and IL-12, leading to the activation of STAT1/STAT4, which in turn activate the T box expressed in T cells (Tbet). Differentiation into Th2 is stimulated by the activation of the GATA binding protein 3 (GATA3) TF through STAT6 signaling. The differentiation of naive T cells into Th17 is governed by the retinoic acid receptor-related orphan receptor gamma t (RORγt) TF, and by the cytokines i.e., IL-6, IL-21, IL-23 and the transforming growth factor beta (TGF-β) (Aggarwal et al., 2003; Harrington et al., 2005; Park et al., 2005; Tesmer et al., 2008). In addition, the TGF-β inhibits T cell differentiation to both the Th1 and Th2 lineages and is also conducive to the cell’s commitment to the iTregs lineage (Schmitt and Williams, 2013).

The complexity of biochemical networks underlying the regulation of T cell differentiation leads to additional questions regarding the mechanisms of the immune response. For instance, based on a large number of possible combinations of extracellular cues, we may ask the following questions: (i) How does the cell decide into which subsequent lineage to differentiate? (ii) What specific combinations of signals are driving a possible switch to a different lineage? (iii) What specific mechanisms are responsible for the T cell differentiation capacity and plasticity?

While regulation of T cell differentiation in the context of the diverse cytokine microenvironment has been studied extensively, effects of the interplay among multiple cytokines on T cell differentiation remain an open question. A systems-level computational model can be used to explore whether, and to what extent, the extracellular cytokine milieu affects the T cell differentiation program. Recently, computational models using various types of mathematical approaches investigated the regulation of phenotypic plasticity, and dynamics in response to diseases (Naldi et al., 2010; Carbo et al., 2013, 2014; Abou-Jaoudé et al., 2014; Martinez-Sanchez et al., 2015). Predictions from these models include novel T cell differentiation pathways (Naldi et al., 2010), transition among T cell types under various micro-environments and perturbations (Martinez-Sanchez et al., 2015), peroxisome proliferator-activated receptor gamma-dependent regulation of Th17 to iTreg switch (Carbo et al., 2013), and IL-21-dependent modulation of IL-10 (Carbo et al., 2014).

Here, we explored the effect of the interplay among extracellular cytokines on differentiation of T cells and their plasticity. We have developed a logic-based computational model (Helikar and Rogers, 2009; Helikar et al., 2012a,b, 2013; Naldi et al., 2015; Abou-Jaoudé et al., 2016; Barberis et al., 2017; Linke et al., 2017) of a signal transduction network that regulates the differentiation process of naive T cells to Th1, Th2, Th17, and iTreg cells and analyzed its dynamics. Local protein–protein regulatory information was manually curated to construct the mechanistic model that contains lineage-specifying TFs (Tbet, GATA3, RORγt, and Foxp3), various signal transducers and activators of transcription (STATs), and other signaling molecules. The model consists of 96 regulatory interactions among 38 components. To explore the entire cytokine microenvironment, we analyzed the model’s dynamics under (i) all possible combinations of extracellular signals, and (ii) various input dosages. The analysis of the model resulted in dynamic signatures that represent previously described, as well as novel cellular phenotypes. These include four canonical phenotypes of differentiated T cells (Th0, Th1, Th2, and iTreg) as well as six complex phenotypes, whereby multiple lineage-specifying TFs are co-expressed. Our results also suggest that the input dosage regulates the balance of specific T cells within the complex T cell phenotypes, providing new insights into specific patterns of environmental input composition and dosage effects on T cell differentiation.

Results

Mechanistic Logical Model of T Cell Differentiation

A comprehensive mechanistic, logic-based model of T cell differentiation was constructed using regulatory information from published literature. The model includes 38 components and 96 biochemical interactions that regulate the differentiation process of major T cell subtypes, such as Th1, Th2, Th17, and iTreg cells. The individual components of the model represent lineage-specifying TFs (Tbet, GATA3, RORγt, and Foxp3), STAT proteins, cytokines, their receptors, and other signaling molecules. The extracellular environment is represented in the model by eight cytokines and a (generic) TCR ligand, known to play a role in T cell differentiation. The network representation of the model is visualized in Figure 1. The regulatory interactions in the model are defined as Boolean functions, which are composed of the “AND,” “OR,” and “NOT” operators (Supplementary Datasheets 1 and 2). The fully annotated model is available for download in a number of formats [including SBML-qual (Chaouiya et al., 2013)], as well as for viewing, and performing simulations, analyses, and additional modifications within the Cell Collective modeling platform¹ (Helikar et al., 2012b, 2013). The model can be accessed directly at: https://www.cellcollective.org/#6678/cd4-t-cell-differentiation.

FIGURE 1

FIGURE 1. Network diagram of the logical model of signal transduction in CD4⁺ T cells. The modeled pathways reflect the canonical pathways known to regulate T cell differentiation into the major effector subtypes (Th1, Th2, Th17) and the regulatory subtype (iTreg). The model includes 38 components, including four lineage-specifying TFs (Tbet, GATA3, RORγt, and Foxp3) and nine extracellular components: TCR-ligand, IFN-γ, TGF-β, IL-4, IL-6, IL-12, IL-18, IL-23, and IL-27. Green arrows represent activation, red arrows represent inhibition, and gray arrows represent the conditions associated with activatory or inhibitory interactions. IFNg_e, interferon-γ (external); IL12, interleukin 12; IL18, interleukin 18; IL23, interleukin 23; IL27, interleukin 27; IL4_e, interleukin 4 (external); IL6_e, interleukin 6 (external); TCR, T cell receptor; TGFb, transforming growth factor-β; Foxp3, Forkhead box P3; GATA3, GATA-binding protein 3; IFNg, interferon-γ; IFNgR, interferon-γ receptor (generic); IL12R, interleukin 12 receptor (generic); IL17, interleukin 17; IL18R, interleukin 18 receptor 1; IL2, interleukin 2; IL21, interleukin 21; IL21R, interleukin 21 receptor; IL23R, interleukin 23 receptor; IL2R, interleukin 2 receptor; IL4, interleukin 4; IL4R, interleukin 4 receptor; IL6, interleukin 6; IL6R, interleukin 6 receptor; IRAK, interleukin-1 receptor associated kinase 1; Jak1, Janus kinase 1; NFAT, nuclear factor of activated T cells 5, tonicity-responsive; NF-κB, nuclear factor of kappa light polypeptide gene enhancer in B cells (generic); RORgt, RAR-related orphan receptor C; SOCS1, suppressor of cytokine signaling 1; STAT1, signal transducer and activator 1; STAT3, signal transducer and activator 3; STAT4, signal transducer and activator 4; STAT5, signal transducer and activator 5; STAT6, signal transducer and activator 6; Tbet, T-box expressed in T cells; TGFbR, transforming growth factor-β receptor (generic).

The model was validated to ensure that it can reproduce differentiation into four canonical phenotypes (Th1, Th2, Th17, and iTreg), as a result of cytokine stimulation and TCR activation (Supplementary Table 1). Furthermore, the model was able to reproduce more complex behaviors (Figure 2). For example, Becskei and Grusby (2007) studied the synergistic effect of positive feedback loops on the expression of the IL-12 receptor (IL-12R). They showed that the number of IFN-γ positive cells and the expression of IL-12R increased when induced by the combination of IL-12 and IL-27. As shown in Figures 2A,B, simulations of the presented model under similar experimental conditions resulted in the same qualitative behavior. Furthermore, it has been experimentally shown that the IL-6 regulates the balance between iTreg and Th17 cells in a dose-dependent manner (Yang et al., 2008; Kimura and Kishimoto, 2010). Similarly, simulations of the model show a clear distinction between iTreg and Th17 in an IL-6-dependent manner (Figure 2C). Finally, simulations of the model, under environmental conditions similar to those that have been shown to induce the mixed Th1–Th2 behavior (Peine et al., 2013), also resulted in a complex phenotype with activation of both Tbet and GATA3 TFs (Figure 2D).

FIGURE 2

FIGURE 2. Validation of the model using complex behaviors. (A–B) Simulations of the model under experimental conditions using the same concentration ratio of IL-27 and IL-12 (i.e., 1000, 100, 10, 1 for IFN-γ and 100, 10, 1 for IL-12R) as in Becskei and Grusby (2007) show that combination of IL-27 with IL-12 leads to a synergistic effect on level of IFN-γ production and the activity of the IL-12 receptor (IL-12R). (C) Simulated IL-6 dose-response effect on the activation of Foxp3 and RORγt. (D) Model simulations reproduce a mixed Th1/Th2 phenotype with varying levels of Tbet and GATA3. The significant activity levels of both Tbet and GATA3 are observed in the area inside the red box.

Novel T Cell Phenotypes Are Predicted by Logical Modeling

With the validated model in hand, we sought to understand its capacity to represent various T cell phenotypes. By using ergodic set analysis [see the section “Materials and Methods” and Todd and Helikar (2012)], we explored the state space of the model under 512 possible combinations of the extracellular stimuli (input compositions) (Figure 3A).

FIGURE 3

FIGURE 3. Simulation of T cell model under all possible environmental conditions. (A) A schematic diagram showing simulation strategy using all possible input compositions (IC). (B) A total of 101 outputs (ergodic sets) were obtained. The number of input compositions stimulating each ergodic set is ranging from lowest 1 to highest 51. Only one ergodic set resulted from a single input composition. (C) Ten phenotypes based on presence and absence of lineage-specifying TFs were obtained. The number of input compositions and ergodic sets for each phenotype are provided in the first and second column. Based on the presence and absence of TFs, each phenotype was determined. Blue cells correspond to TFs present in the phenotype, whereas the orange cells represent inactive TFs.

A total of 101 ergodic sets (outputs) were obtained as a result of 511 input compositions (Supplementary Table 2). Out of the 511 compositions, 45 input compositions resulted into fixed points (a single remaining input composition was not analyzable even on a supercomputer due to the large size of state space which could not be computed on a feasible temporal scale). The number of input compositions for each output ranged from 1 to 51 (Figure 3B). We obtained one output (output 3) that can be stimulated by the maximum of 51 input compositions. Two outputs (outputs 6 and 13) were each stimulated by the maximum of 48 input compositions (Figure 3B). Furthermore, four outputs (outputs 10, 22, 29, and 32) were each achieved by 16 different input compositions. All outputs that are individually stimulated by 16 or more input compositions have input compositions with an inactive TCR ligand.

The number of input compositions for the remaining outputs varied from 1 to 4. These input compositions contained an active TCR ligand. In this group, a total of 37 outputs were obtained, whereby each of them was stimulated by four input compositions. A total of 56 outputs were each stimulated with two input compositions. Only one output was stimulated by a single input composition.

Thus, 7 (out of 101) outputs were achieved when stimulated by 211 input compositions with the absence of a TCR ligand. On the other hand, 94 outputs (out of 101) were obtained when stimulated by 255 input compositions with an active TCR ligand. Therefore, fewer outputs (101) have been observed than the total number of input compositions (511), suggesting that a specific cell fate (output) can result from multiple signal compositions, processed by a cell based on biochemical rules in a signaling network (Helikar et al., 2008; Balázsi et al., 2011; Palau-Ortin et al., 2015).

Next, we explored the biological relevance of the produced outputs. As the model centers on the regulation of T cell phenotypes and the TFs related to each differentiated T cell subtype, we classified all the outputs based on the presence of the four TFs (GATA3, Tbet, RORγt, and Foxp3). We found that the model outputs (as a result of the 511 input compositions) cluster into 10 biologically relevant phenotypes. These include the canonical (single cell type) phenotypes as well as the complex phenotypes having more than one lineage-specifying TF.

Specifically, we found four canonical T cell phenotypes that carried Tbet, GATA3, or Foxp3, representing Th1, Th2, and iTreg, respectively (Figure 3C). Furthermore, we found that 219 input compositions resulted in nine outputs with no TFs present (Th0 phenotype). Most of the outputs that represent the Th0 phenotype (>95%) were stimulated by the input compositions with an inactive TCR ligand. The remaining Th0-leading input compositions contained an active TCR ligand along with IL-23, or IL-18, or IL-6. This corresponds to the experimentally established scenarios, whereby the T cells cannot differentiate in the absence of TCR activation or in the absence of key lineage-specific cytokines (Podojil and Miller, 2009; Zhu et al., 2010; Chen and Flies, 2013). Fifty-two input compositions led to 16 outputs with active Tbet, representing the Th1 phenotype. A total of 24 input compositions produce 10 outputs with active GATA3, representing the Th2 phenotype, while four input compositions led to one output with active Foxp3, representing the iTreg phenotype. We did not observe distinct outputs with only RORγt active; instead, RORγt was part of the complex phenotypes (discussed below).

In addition to the four canonical phenotypes, the model predicted six complex phenotypes. The number of input compositions, and the number of outputs that represent each complex phenotype, is summarized in Figure 3C. Of the six complex phenotypes, three of them including Th1–Th2 (Hegazy et al., 2010; Evans and Jenner, 2013; Peine et al., 2013), Th1–iTreg (Koch et al., 2009), and Th17–iTreg (Eisenstein and Williams, 2009) were experimentally observed earlier, thus further validating the model. The model also predicted three novel complex phenotypes, Th1–Th2–iTreg, Th1–Th17–iTreg, and Th1–Th2–Th17–iTreg, for which experimental validation is foreseeable.

Cytokine Composition Establishes T Cell Phenotypes

Once the classification of all the model outputs into biologically relevant phenotypes was carried out, we analyzed the input compositions (environmental conditions) leading to each of the 10 biological phenotypes. This analysis resulted in 27 patterns of input compositions (Figure 4). We also identified the minimal input compositions that are needed to stimulate each phenotype (Figure 5). Additionally, the signal transduction sub-networks activated for each phenotype, simulated under a representative input composition, are shown in Figure 6.

FIGURE 4

FIGURE 4. Input compositions for all T cell phenotypes. The color map shows the patterns of input compositions that give rise to each observed phenotype. For example, co-expression of RORγt–Foxp3 needs TCR + IL-6 + TGF-β (+IL-18 and IL-23 can also be active) to be active and IL-12, IFN-γ, IL-27, and IL-4 to be inactive.

FIGURE 5

FIGURE 5. Minimal input compositions required to stimulate T cell phenotypes. All inputs in blue boxes are required to stimulate the corresponding phenotype. For example, the Th1–Th2 phenotype can be stimulated by input composition: TCR ligand + (IFN-γ OR IL-12 OR IL-27) + (IL-4 OR IL-12 + IL-18 + IL-6). Minimum three inputs are required to stimulate the Th1–Th2 phenotype (e.g., TCR ligand + IFN-γ/IL-12/IL-27 + IL-4).

FIGURE 6

FIGURE 6. Active signal transduction sub-networks for T cell phenotypes. The active sub-networks under an input composition are mapped on the entire network of T cell differentiation and are marked by bold arrows. Green arrows represent activation and red arrows represents inhibition, and gray arrows represent the condition associated with activation and inhibition. (A) Th1 phenotype under input composition: TCR-ligand + IL-12 + IL-27 + IFN-γ, (B) Th2 phenotype under input composition: TCR-ligand + IL-4, (C) iTreg phenotype under input composition: TCR + TGF-β, (D) Th17-iTreg phenotype under input composition: TCR-ligand + IL-6 + TGF-β, (E) Th1–Th2 phenotype under input composition: TCR-ligand + IL-4 + IL-12 + IL-27 + IFN-γ, (F) Th1–iTreg phenotype stimulated under input composition: TCR-ligand + IL-12 + IL-27 + IFN-γ + TGF-β, (G) Th1–Th17–iTreg phenotype stimulated under input composition: TCR-ligand + IL-12 + IL-27 + IFN-γ + IL-6 + TGF-β, (H) Th1–Th2–iTreg phenotype stimulated under input composition: TCR-ligand + IL-4 + IL-12 + IL-27 + IFN-γ + TGF-β, (I) Th1–Th2–Th17–iTreg phenotype stimulated under input composition: TCR-ligand + IL-4 + IFN-γ + IL-6 + TGF-β.

As indicated in the model validation section, we found that the canonical phenotypes (Th0, Th1, Th2, and iTreg) are regulated by one or more cytokines. We also found that all the complex phenotypes can be stimulated by more than one input composition. Strikingly, our modeling effort predicts that in order to induce specific phenotypes, certain cytokines cannot be co-present in a given input composition (Figures 4, 5). For example, based on our model, TGF-β should not be present in the input compositions leading to the Th1–Th2 phenotypes, and IL-6 should be absent from the input compositions that lead to iTreg, Th1–iTreg, and Th1 phenotypes. On the other hand, IL-4 can be present in the input composition leading to Th1, but only when co-present with IL-6. IL-4 also needs to be absent in input compositions leading to iTreg, Th17–iTreg, Th1–iTreg, and Th1–Th17–iTreg phenotypes. Finally, IL-12 and IL-18 cannot be co-present in the extracellular environment that stimulates differentiation into Th1, Th2, Th1–iTreg, Th1–Th17–iTreg, and Th0 (in the absence of the TCR ligand) phenotypes.

The previously mentioned heterogeneous and conditional effect of combining IL-12 and IL-18 is also supported and partially explained through experimentally described regulatory mechanisms (Yoshimoto et al., 1998; Nakanishi et al., 2001). Specifically, we observed that combining IL-18 with IL-12 favors co-expression of Tbet, GATA3, and Foxp3. It was previously shown that combining IL-12 and IL-18 can synergistically increase the Tbet-stimulated IFN-γ production in Th1 cells (Tominaga et al., 2000). In another study, it was shown that IL-18, but not IL-12, increases the production of IFN-γ by CD8+ and CD4⁺ T cells in the K14E7 transgenic skin (Gosmann et al., 2014). Further, the combination of IL-12 and IL-18 has been shown to induce the production of IFN-γ in the absence of antigen (Munk et al., 2011). Finally, it has been shown that IL-18 in the absence of IL-12 can stimulate Th2 response (Nakanishi et al., 2001).

To further investigate the effect of IL-12 and IL-18 on the Th1–Th2–iTreg phenotype, the model was simulated under the input composition of IL-12, IL-18, and TCR (with all other cytokines inactive). Simulation results suggested the synergistic effect of IL-12 and IL-18 on the activity level of GATA3 and Foxp3. Interestingly, the increased activity of GATA3 and Foxp3 was observed in the absence of external IL-4 and TGF-β (Figure 7A), suggesting that the combination of IL-12 and IL-18 (while controlling for the TCR signal) are able to stimulate the Th1–Th2–iTreg phenotype in an IL-4- and TGF-β-independent manner. We also found that the combination of IL-12 and IL-18 is a weaker activator of GATA3 and Foxp3 (Figure 7A). This is because the IL-12 can also stimulate Tbet, which in turn suppresses the GATA3 and Foxp3. Results obtained from the simulated IL-12R knock-out suggested an eightfold increase in the activity of GATA3, whereas the overexpression of IL-12R slightly decreased the activity levels of GATA3 and Tbet. Knock-out of IL-18R resulted in a complete inactivation of GATA3 and Foxp3, whereas the overexpression of IL-18R resulted in a greater than twofold increase in the activity levels of GATA3 and Foxp3. These results indicate that the knock-out of IL-12R favors Th2 phenotype, whereas the knock-out of IL-18R favors Th1 phenotype under Th1–Th2–iTreg stimulating environmental conditions.

FIGURE 7

FIGURE 7. Combination of IL-12 and IL-18 favors the Th1–Th2–iTreg phenotype. (A) The model was simulated using input composition (TCR ligand + IL-12 + IL-18) under five conditions, including wild-type, and knock-out and overexpression of both IL-12R and IL-18R. Overexpression of IL-12R favors Th1 phenotype, whereas IL-12R knock-out significantly increased the activity level of GATA3. IL-18R knock-out completely inactivated both GATA3 and Foxp3. (B) A subnetwork that might induce GATA3 and Foxp3 in the presence of IL-12 and IL-18, but in the absence of IL-4 and TGF-β. (C) Knock-out of STAT5 and NF-κB resulted in complete inactivation of GATA3 and Foxp3, whereas overexpression of STAT5 induced a strong Foxp3 response. Overexpression of NF-κB increased the activity of both GATA3 and Foxp3 by more than twofold.

The differentiation to Th2 was previously observed in airway epithelia in the presence of IL-18, but not IL-4 (Murai et al., 2012). The IL-4-independent Th2 stimulation possibly occurs through the STAT5-mediated GATA3 activation (Yamane et al., 2005; Paul, 2010). The IL-18R1 signaling was also found to promote Foxp3+ iTreg cell function within colonic lamina propria (Harrison et al., 2015). To better understand the mechanism of how the IL-18 and IL-12 can stimulate GATA3 and Foxp3, we further analyzed the network structure of the model. We found that IL-12 and IL-18 can possibly induce the production of IL-2, which stimulates GATA3 and Foxp3 in STAT5-dependent pathways (Figure 7B). The knock-out simulation of NF-κB or STAT5 resulted in complete inactivation of GATA3 and Foxp3. On the other hand, the overexpression of STAT5 increased the mean activity level of Foxp3 by 62-fold, while no change in activity levels of GATA3 was observed. The simulated over-expression of NF-κB had shown 5.4-fold and twofold increase in the activity levels of Foxp3 and GATA3, respectively. These results predict the role of IL-12 and IL-18 in stimulation of the Th1–Th2–iTreg phenotype in an NF-κB- and STAT5-dependent manner (Figure 7C). Furthermore, our simulation results suggest that a combination of IL-18 and IL-12 can stimulate Tbet, GATA3, and Foxp3; however, the activity levels of GATA3 and Foxp3 were lower than that of Tbet (Figure 7A). Additionally, we have found that IL-12 and IL-18 combination in the presence of IL-6 can stimulate the Th1–Th2 phenotype (Supplementary Figure 1 in Supplementary Datasheet 3).

Altogether, we have identified input composition patterns that include the minimum combinations of cytokines required to stimulate a particular T cell phenotype, as well as complete pattern of cytokines that can be co-present to stimulate a given phenotype (Supplementary Table 3). Our results also predict the relevance of IL-12 and IL-18 in regulating the Th1–Th2–iTreg phenotype. Finally, we predicted an alternative pathway that can stimulate GATA3 and Foxp3 in an IL-4 and TGF-β-independent manner.

Cytokine Dosage Determines the Balance Between Complex T Cell Phenotypes

In the previous section, various input compositions that lead to different canonical and complex phenotypes were characterized. The logical question that we raise now is: How is the balance of each T cell subtype within a complex phenotype controlled?

As indicated in the “Introduction” section, several reports suggest that the balance between Th17 and iTreg is regulated by the dosage of IL-6 (Kimura and Kishimoto, 2010; Omenetti and Pizarro, 2015). To explore how the input dosages within each composition affect the complex phenotypes, we analyzed the model under various activity levels of cytokines and the TCR ligand under the complete set of input compositions.

We used the representative input compositions for each identified phenotype as described in Figures 4, 5. Specifically, we used two types of representative input compositions from each row in Figure 4. The two types include, one with the maximum number of inputs that can be simultaneously present to stimulate a specific T cell phenotype, and a second type that is represented by input compositions consisting of the minimal number of inputs required to stimulate the identified phenotypes.

Sensitivity analysis of the model was performed to describe the effect of each input in its composition on the TF(s) for the corresponding complex phenotype (Figure 4). The analysis predicted individual inputs that are important for regulating the balance among lineage-specifying TFs. For example, for the Th1–Th2 phenotype, when stimulated with a maximum of eight inputs, the sensitivity analysis suggested that IL-27, IFN-γ, and IL-12 are negatively correlated with GATA3 (Figures 8A–C). The TCR signal is negatively correlated with GATA3 [partial correlation coefficient (PCC) range = -0.18 to -0.19] under three input compositions (Figures 8A–C). Interestingly, a positive correlation between the TCR ligand and GATA3 was observed when the Th1–Th2 phenotype was stimulated in the absence of IL-4 (and in the presence of IL-12, IL-18, and IL-6) (Figure 8D). On the other hand, the IL-18 had a moderate negative correlation with Tbet (PCC range = -0.28 to -0.29) under all tested input compositions (Figures 8A–D). The IL-4 had a very low correlation with Tbet (PCC range = 0.005–0.01) under all tested input compositions (Figures 8A–D).

FIGURE 8

FIGURE 8. Sensitivity analysis showing the input effect on lineage-specifying TFs for the Th1–Th2 phenotype. Panels (A–D) are based on simulations using maximal input compositions. Panels (E–H) are based on minimal input compositions. PCC as a measure of association between inputs (cytokines and TCR) and lineage-specifying TFs is shown on Y-axes and input composition (cytokines and TCR) is shown on X-axes.

Next, in the case of the Th1–Th2 phenotype stimulated under minimal input compositions, higher correlations between the inputs and TFs were observed compared to the correlations observed with maximal input compositions (Figures 8E–H). In the case of the input composition “TCR + IL-12 + IL-4,” IL-12 had a strong negative correlation (PCC = -0.68) with GATA3, and a strong positive correlation with Tbet (PCC = 0.7). In this case, the TCR ligand had a moderate negative and positive correlation with GATA3 (PCC = -0.25) and Tbet (PCC = 0.22), respectively. In the case when Th1–Th2 phenotype was stimulated under input composition “TCR + IL-12 + IL-18 + IL-6,” the TCR ligand was positively correlated with both GATA3 (PCC = 0.30) and Tbet (PCC = 0.30). A strong positive correlation was observed between IL-4 and GATA3 (PCC = 0.65) under the input composition “TCR + IFN-γ + IL-4.” In the Th1–Th2 complex phenotype, we observed that the TCR ligand is negatively correlated with GATA3. The negative effect of a strong TCR ligand signal on GATA3 is in agreement with the earlier studies suggesting that a strong TCR signal can promote a strong Th1 response, whereas a weaker signal favors the Th2 response (van Panhuys et al., 2014). The sensitivity analysis results for all other mixed phenotypes are provided in Supplementary Table 4.

In summary, the sensitivity analysis of our model predicts “driver” inputs. Furthermore, it characterizes the strength and direction (positive or negative) of the effect inputs can have on the regulation of the balance of each T cell subtype within the complex phenotypes. The strength of association between the inputs and TFs varied based on the number of inputs in the input compositions.

Determining the Optimal Input Dosage Regulating the Balance Between Complex Phenotypes

In the previous sections, the predicted complex T cell phenotypes, input compositions, as well as the potential dosage effect each input can have on the phenotype, were discussed. Next, we examined the specific activity levels of the input compositions required to control each specific T cell phenotype. The model was simulated under 10,000 randomly generated environmental conditions within the context of each relevant input composition. Results from these simulations provided us with specific input activity levels that have a low coefficient of variance (CV) in activity levels of co-expressed lineage-specifying TFs. Specifically, we investigated and characterized the activity levels for each input composition that will drive a complex T cell phenotype to each of the T cell subtypes or a balanced mixed phenotype by maximizing the activity levels of the respective TFs. For example, to achieve a balanced Th1–Th2 phenotype that has similar activity levels to that of GATA3 and Tbet, we characterized the optimal activity levels for each input in the Th1–Th2-leading input compositions. This predicted optimal input composition includes the low activity of the TCR ligand, IFN-γ, IL-12, and IL-27, medium activity of IL-18 and IL-6, and high activity of IL-4. The activity of IL-23 can vary from low to high, whereas TGF-β should be inactive (Figure 9A).

FIGURE 9

FIGURE 9. Dose-dependent regulation of complex phenotypes co-expressing lineage-specifying TFs. (A) Optimal settings using maximal input composition that stimulate co-expression of lineage-specifying TFs. Shown are median activity levels of inputs that result in a balanced activity level of both TFs. Activity levels ranging from 0 to 40 (blue) were considered as low, 40 to 80 were medium (light blue), and 80 to 100 were high (purple). (B) Activity levels of GATA3 and Tbet when simulated under optimal input settings that stimulate the Th1–Th2 phenotype. Median activity levels from panel A were used for the simulation. (C–J) Dose-response curves illustrating the effect of TCR, IL-12, IL-18, and IL-27 on GATA3 and Tbet.

To illustrate the effect of using optimal activity levels, we stimulated the Th1–Th2 phenotype by using a median value of optimal activity level for each input. As expected, the simulation results show similar activity levels of Tbet and GATA3 (Figure 9B). To further investigate the effect of dominant inputs (identified from the sensitivity analysis) on the Tbet–GATA3 combination, we performed dose-response analysis by varying the dominant cytokines while fixing the other inputs to median activity levels (Figure 9A). As expected, our results (Figures 9C–J) suggest that the increased signal strength of TCR ligand or increased activity of IL-12 and IL-27 can drive the Th1–Th2 phenotype toward Th1 by increasing the activity of Tbet and decreasing the activity of GATA3. In contrast, the increased activity of IL-18 can drive Th1–Th2 phenotypes toward Th2.

Discussion

In this study, we sought to investigate the cellular phenotypes as a result of CD4⁺ T cell differentiation under diverse environmental conditions and understand how the balance between complex phenotypes is regulated. To achieve this, by manually curating literature data, we constructed a mechanistic computational (logical) model of signal transduction that regulates the differentiation of naive T cells into Th1, Th2, Th17, and iTreg cells. The components (i.e., proteins and genes) in a logical model can have binary (0 or 1) states at any time t. The state of the network evolves stepwise based on the logical rules defined for each model component (Helikar and Rogers, 2009; Helikar et al., 2012a,b, 2013; Naldi et al., 2015; Abou-Jaoudé et al., 2016; Barberis and Verbruggen, 2017; Linke et al., 2017).

We systematically characterized the model’s dynamics in the context of activity of lineage-specifying TFs under 511 input compositions consisting of eight cytokines and a TCR signal. In addition to the dynamics representing the classical Th0, Th1, and Th2 phenotypes, we found several complex (mixed) phenotypes (dynamics with more than one lineage-specific TFs), including Th1–Th2, Th1–iTreg, Th17–iTreg, Th1–Th2–iTreg, Th1–Th17–iTreg, and Th1–Th2–Th17–iTreg. Our results are in agreement with recent studies that reported hybrid T cell phenotypes in vitro and in vivo (Zhou et al., 2008; Peine et al., 2013). Stable complex Th1–Th2 phenotypes parallel to the classical Th2 phenotypes were observed in vivo upon infection mediated by parasites Schistosoma mansoni and Heligmosomoides polygyrus (Peine et al., 2013), as well as by the threadworm Strongyloides stercoralis (Bock et al., 2017). Moreover, Th1–iTreg intermediate phenotypes were observed during Th1 polarizing infections (Koch et al., 2009; Oldenhove et al., 2009; Evans and Jenner, 2013). In a recent system level study, a continuum of T cell differentiation states with stable co-expressed lineage-specific TFs has been observed when stimulated under different combinations of six cytokines (Eizenberg-Magar et al., 2017).

Interestingly, we did not observe a canonical Th17 (RORγt-only) phenotype. Instead, our model predicts the existence of a mixed Th17–iTreg phenotype. This result can be partially explained by the fact that both Th17 and iTreg share a common mechanism by cytokine TGF-β, and the differentiation of naive T cells into iTreg or Th17 depends on the cytokine-driven (TGF-β and IL-6) balance of lineage-specifying TFs Foxp3 and RORγt (Omenetti and Pizarro, 2015). In addition, it is known that the Th17/Treg balance is critical to maintain immune tolerance. The imbalance of Th17/Treg has been observed in the peripheral blood of cervical cancer patients (Chen et al., 2013), non-small cell lung cancer patients (Duan et al., 2015), and in patients with chronic low back pain (Luchting et al., 2014). Thus, the complex Th17–iTreg phenotype might play an important role in maintaining Th17/Treg homeostasis. Such complex RORγt–Foxp3 co-expressing T cells were observed in an autoimmune diabetes model (Ichiyama et al., 2008; Tartar et al., 2010), in the lamina propria (Zhou et al., 2008), in the peripheral blood and tonsils (Voo et al., 2009), and in the large intestine (Ohnmacht et al., 2015; Fang and Zhu, 2017). It is also possible that the lack of Th17-only phenotype is due to the incomplete nature of the model. However, it suggests that additional experimental validation may be required to better understand the relationship and mechanism of switching between iTreg and Th17 phenotypes.

We also predicted novel phenotypes that have the potential to have three active TFs (Tbet–GATA3–Foxp3, Tbet–RORγt–Foxp3), as well as one with all four TFs (Tbet–GATA3–RORγt–Foxp3). In partial support of our prediction, basal levels of Tbet and GATA3 have been observed in iTreg cells (Yu et al., 2015). While not yet shown experimentally, the Th1–Th17–iTreg phenotype was also predicted by a similar modeling approach (Naldi et al., 2010).

By analyzing all possible inputs combinations, we obtained the minimal and maximal input compositions for each identified phenotype. The minimal composition includes a minimum number of inputs that can stimulate a phenotype. On the other hand, the maximal composition includes a maximum number of inputs that can be simultaneously active to result in the same phenotype. In this analysis, we found that in order to stimulate Th1, Th2, Th1–iTreg, Th1–iTreg, Th1–Th17–iTreg, and Th0 phenotypes, IL-12 and IL-18 cannot be combined in the environment. We observed that the combination of IL-12 and IL-18 leads to the stimulation of GATA3 and Foxp3 even in the absence of IL-4 and TGF-β via a NF-κB-dependent pathway. We predicted that a combination of IL-18 and IL-12 could result in a Th1–Th2–iTreg complex phenotype. Analysis of the model’s network structure suggests a potential mechanism that is dependent on NF-κB and STAT5 (Figure 7B). Previous studies suggest that IL-18 has a context-specific functional heterogeneity and can induce both Th1 and Th2 T cell phenotypes. The combination of IL-12 and IL-18 has been shown to have a synergistic effect on IFN-γ production that stimulates the Th1 phenotype (Tominaga et al., 2000; Munk et al., 2011). It has also been found that IL-18 alone (without IL-12) can stimulate the Th2 phenotype (Nakanishi et al., 2001). In a study on airway epithelial cells in response to Alternaria, it was found that secreted IL-18 has the capacity to stimulate the Th2 phenotype (Murai et al., 2012). Since IL-12 can up-regulate IL-18R expression, it might be possible that the combination of IL-12 and IL-18 may regulate the Th1, Th2, Th1–Th2, and Th1–Th2–iTreg phenotypes in a dose-dependent manner.

Next, the sensitivity analysis of the model suggested that the dosage of the individual inputs regulates the balance within the different complex T cell phenotypes. We investigated the dosage effect by using both minimum and maximum number of inputs under varying activity levels. For example, our results suggest that the dynamics of the complex Th1–Th2 phenotype depend on the combination and dosage of IFN-γ, IL-12, IL-27, IL-18, IL-4, and the TCR ligand. The increased activity levels of the cytokines IFN-γ, IL-12, IL-27, and TCR ligand drive the phenotype toward Th1, whereas the IL-18 or IL-4 drive the Th2 phenotype. The IL-23 and IL-6 have no correlation with either Tbet or GATA3. Under both maximal and minimal input compositions, the IL-4 had low to no correlation with Tbet. On the other hand, the IL-18 was positively correlated with GATA3 and negatively correlated with Tbet. Thus, we predicted that IL-18 may have a dominant role over IL-4 to favor Th2 phenotype under the Th1–Th2 stimulating environmental conditions.

Next, we identified the activity levels of the inputs required to regulate the complex T cell phenotypes. Our results suggest a range of activity levels required to obtain a specific phenotype under minimal and maximal input compositions. For example, a high amount of IL-4 or IL-18 and a low amount of IFN-γ, IL-12, IL-27, and TCR ligand are required to stimulate the Th1–Th2 phenotype under maximal input composition. Low activity of GATA3 under higher TCR ligand activity is indeed in agreement with the literature where it has been shown that a strong TCR signal represses GATA3 (Aguado et al., 2002; Yamane et al., 2005; Paul, 2010; Altin et al., 2011; Yamane and Paul, 2012). Interestingly, our results showed an increase as well as a decrease in the activity levels of GATA3 depending on the activity levels of IL-12. This can be achieved as a result of IL-12 up-regulating IL-18R, which induces NF-κB-mediated GATA3 activation. On the other hand, a higher activity of IL-12 results in a strong Tbet activation, which in turn suppresses GATA3. Although the predicted activity levels are dimensionless and semi-quantitative, they provide a starting point for calibrations against ligand concentrations in specific experimental research protocols.

In summary, results provided in this study can provide a platform to generate and design testable hypotheses in the context of T cell differentiation in response to various combinations and dosage of environmental signals. Furthermore, the presented results and the mechanistic model can be used as tools to further investigate the specific pathway mechanisms that govern each complex phenotype. Input availability and relative dosage at which inputs generate a productive signaling cascade necessarily result in a variable timing of an immune response. Specifically, we and others propose that dosage- and timing-dependent impact of inputs, such as ILs, may impact the T cell differentiation (Barberis et al., 2018; Martinez-Sanchez et al., 2018). This may be investigated by employing experimental methodologies that we have recently envisioned (Barberis and Verbruggen, 2017). Furthermore, crosstalk between ILs and signaling cascades, such as the one governing the cell cycle, may impinge on a timely T cell-mediated protective response (Barberis et al., 2018). These aspects are the focus of our current research efforts. Together with new model-based predictions, improving the understanding of the detailed mechanisms underlying T cell differentiation, can be helpful to design strategies for immunotherapy against pathogens and various diseases of the immune system.

Materials and Methods

Model Construction

The computational model is a mechanistic, logic-based model of signal transduction processes known to regulate CD4⁺ T cell differentiation into Th1, Th2, Th17, and iTreg cells. Each component of the model can assume an active (1) or inactive (0) state at any time t. The activity state of the model’s internal components is determined by the regulatory mechanisms of other directly interacting components. These regulatory mechanisms are described with Boolean functions (Samaga and Klamt, 2013; Albert and Thakar, 2014; Le Novère, 2015; Naldi et al., 2015; Abou-Jaoudé et al., 2016; Linke et al., 2017).

The new signal transduction model was constructed manually by curating published regulatory mechanisms of each signal transduction component. Each of the 38 components in the model corresponds to a signaling molecule (mainly proteins). The model also contains nine external components that represent the extracellular environment, consisting of eight cytokines (IFN-γ, TGF-β, IL-4, IL-6, IL-12, IL-18, IL-23, and IL-27) and a generic TCR ligand. The final model consists of 38 components (29 internal and 9 external) connected with 96 interactions. The model is fully annotated with published evidence for each component and interaction to ensure transparency and reproducibility. The model is available via the web-based modeling and analysis platform Cell Collective (Helikar et al., 2012b, 2013), accessible at https://www.cellcollective.org (under Published Models) where it can be simulated as well as downloaded (and other logical models published by the community) in several file formats (such as SBML-qual, text file of logical functions, and truth tables).

State Space Analysis

The logical model herein is a Probabilistic Boolean Control Network (PBCN) (Todd and Helikar, 2012), whereby each external input (components that are not regulated by other model components) is activated by a user-defined probability of activation (ranging from 0 to 100%). The activity levels of the external inputs and the logical rules associated with each internal node allow the system to update stochastically in time. As such, a PBCN is a reducible Markov chain (Tijms, 2003). We used ergodic sets (recurrent communicating classes of the corresponding Markov chain) as a model of stable cell states that represent the phenotype of a differentiated T cell. Ergodic sets are a collection of states in state space such that once the system evolves to one of these states it will remain in this set of states. In this way, the ergodic sets are the stochastic equivalents to attractors in purely Boolean networks (Ribeiro and Kauffman, 2007).

From each initial condition, the system will arrive in one of a (possibly) different collection of ergodic sets. In order to find all the ergodic sets, one would need to let the system evolve from every possible initial condition. Given the large number of possible initial conditions (2²⁹), this is computationally infeasible. Thus, we found those ergodic sets that can be reached from the initial state where all internal components are inactive. This represents our goal, i.e., to identify cell phenotypes that are the result of differentiation from naive T cells (i.e., all model components are inactive). Once an ergodic set was identified it was treated as an irreducible Markov chain and thus has an associated limiting distribution. Activities of the internal components are interpreted by approximating the limiting distribution of the Markov chain via simulations in Cell Collective. This means that each internal component has a unit less activity level corresponding to the probability that it is active in the limiting distribution of the Markov chain.

Identification of Ergodic Sets

The extracellular environment (external input components) in the presented model consists of nine stimuli — eight cytokines and a generic TCR ligand. A given extracellular environment is described according to those stimuli that are off (no activity) and those that are on (some level of activity). Thus, there are 2⁹ = 512 possible off/on configurations for the extracellular environment (input compositions). The ergodic sets that are reachable from the naive state (where all components are inactive) depend only on this off/on description and not on the activity level of the non-off cytokines. We were able to identify the corresponding reachable ergodic sets for 508 of these input compositions. The only extracellular environments that are yet unknown are the three where all stimuli are on except for TGF-β, or IL-23, or IL-4. The ergodic sets were identified in two steps. In the first step, Tarjan’s algorithm (Tarjan, 1972) was used to identify communicating classes of states. In the second step, these classes were directly tested to determine if they were closed. The ergodic sets (other than the fixed points) ranged in size from the smallest, with two states, to the largest with 594,962 states. These ergodic sets correspond to the “outputs” in Figures 3A,B. Each state in an ergodic set specifies the state of the internal network. In order to classify an ergodic set, for each internal component we computed the percentage of states in which the component was active. For example, the ergodic set that was identified when the TCR ligand and IL-4 are off while all other external stimuli are on, was found to have 64 states. Each of IL-18R, IL-4R, IRAK, NF-κB, and STAT6 are on in 50% of states, though not the same 50% of states. All other internal components were off in all of the 64 states. In this case, as no lineage-specific TFs are expressed at any level, it is classified as a Th0 phenotype.

The computations to find the ergodic sets were implemented in PERL and were run on an 82-node Linux cluster. Most computations of the ergodic sets required around 10–20 gigabytes of RAM and took from hours to days for the Tarjan’s algorithm to find an ergodic set. (Some required much more). In general, given an initial condition and off/on input composition, several ergodic sets could be reached. We found that out of the 512 possible input compositions, 502 compositions lead to a unique ergodic set and 6 of them lead to exactly two ergodic sets. There were three input compositions that led to one ergodic set, but for which the algorithm had not finished the complete search even after 7 days of calculations. Thus, for these three input compositions, there could be reachable ergodic sets that we did not identify. One input composition, in which all external inputs are active, ran for 7 days without finding any ergodic sets (this is the only input composition for which we have no ergodic set). As we got inconclusive results from the aforementioned incomplete analyses, the corresponding four input compositions were excluded from any reported results.

Model Simulations in Cell Collective

Model simulations were performed in the web-based modeling platform, Cell Collective². Although the model is built by using discrete mathematics, the output activity levels of individual components can be represented as semi-continuous values ranging from 0 to 100% as previously described in Helikar et al. (2008) and Helikar and Rogers (2009). Each simulation was conducted using synchronous updates, and consisted of 5,000 steps, where the activity level of the measured output component was calculated as the fraction of ones (active states) over the last 500 iterations that describe the model’s steady behavior (Helikar et al., 2008; Helikar and Rogers, 2009). The activity levels (dosage) of external components is unit-less and defined as a per-cent chance (probability ^∗ 100) of the component being active in a given time t. Depending on the desired experiment, the activity levels of external components can be set by the user to specific values, or they can be set to ranges from which values during each simulation are selected randomly (e.g., to simulate dose-response experiments).

Once the ergodic sets were identified, expressions of the internal components and their dependencies on the dosages of the external cytokines and the TCR ligand were investigated via the Cell Collective (Helikar et al., 2012b).

For each ergodic set, we chose one of its states as an initial condition and then simulated the model with the corresponding extracellular conditions via the Cell Collective. For each of the active input cytokines, the activity levels varied between 1 and 99%. Further details of the use of the Cell Collective are specific to the types of analysis as described below.

Sensitivity Analysis

The model was simulated in Cell Collective, whereby the activity levels of the inputs for each composition varied. By using the model-generated simulation data under 10,000 randomly generated environmental conditions, the association between inputs (cytokines and TCR ligand) and outputs (lineage-specifying TFs) was determined by probabilistic global sensitivity analysis based on PCC using the “sensitivity” package in R (R Development Core Team, 2011; Pujol et al., 2017). The PCC measures the strength of association between the output and input parameters after removing the linear effect of other input parameters (Marino et al., 2008; Pujol et al., 2017). The PCC between input and output is the correlation coefficients between residuals (x_j - _j) and (y -ŷ), where x_j and y are input and output, respectively, and _j ŷ are linear regression models [shown in Equation (1)] (Marino et al., 2008).

$${\widehat{x}}_j=c_0+{\sum }_{\begin{array}{l}p = 1\hfill \\ p \ne j\hfill \end{array}}^k{c}_p{x}_p\text{and}\widehat{y}=b_0+{\sum }_{\begin{array}{l}p=1\hfill \\ p \ne j\hfill \end{array}}^k{b}_p{x}_p.\left(1\right)$$

Optimal Settings Analysis

Once again, the model was simulated using 10,000 randomly generated environmental conditions for each input composition that can stimulate a complex phenotype. We sought to identify the environmental conditions wherein multiple lineage-specifying TFs can have balanced activity levels. First, we used the CV [Equation (2)] between TFs to measure variability. Further, we selected simulation results under which the lowest variability between TFs was observed. We selected corresponding environmental conditions that had lowest CV among TFs. Next, we selected the top 10 environmental conditions based on the outputs that have the highest activity levels of TFs. Thus, we considered both the balance of activity levels as wells as the quantity of co-expressed TFs. Finally, we defined ranges of activity levels of inputs from the selected environmental conditions. Further, for Th1–Th2, we simulated the effect of dominant inputs by individually varying IL-12, IL-18, IL-27, and the TCR ligand and using median activity levels from identified optimal activity levels for other inputs. We used R-scripts to determine the optimal activity levels from simulation data obtained via Cell Collective (Helikar et al., 2012b). The effect of dominant inputs on TFs in a complex phenotype was shown using the Generalized Additive Model (GAM) fitted scatter plots generated using “ggplot2” package in R.

$$\%\text{CV}\text{=}\frac{\text{Standard}\text{deviation}}{\text{mean}}\times \text{100}.\left(2\right)$$

Author Contributions

BP, RT, DB, MB, and TH designed the research. BP, RT, AM, and DB performed the research. BP, MB, and TH analyzed the data. BP, RT, MB, and TH wrote the paper.

Funding

This project was supported by NIH grant no. 5R35GM119770-02 to TH and by the SILS Starting Grant of the University of Amsterdam to MB.

Conflict of Interest Statement

TH has served as a scientific advisor and/or consultant to Discovery Collective.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We would like to thank Resa Helikar and Robert Moore for providing feedback on the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2018.00878/full#supplementary-material

Footnotes

1. ^ https://www.cellcollective.org
2. ^ https://www.cellcollective.org

References

Abou-Jaoudé, W., Monteiro, P. T., Naldi, A., Grandclaudon, M., Soumelis, V., Chaouiya, C., et al. (2014). Model checking to assess T-helper cell plasticity. Front. Bioeng. Biotechnol. 2:86. doi: 10.3389/fbioe.2014.00086

PubMed Abstract | CrossRef Full Text | Google Scholar

Abou-Jaoudé, W., Traynard, P., Monteiro, P. T., Saez-Rodriguez, J., Helikar, T., Thieffry, D., et al. (2016). Logical modeling and dynamical analysis of cellular networks. Front. Genet. 7:94. doi: 10.3389/fgene.2016.00094

CrossRef Full Text | Google Scholar

Aggarwal, S., Ghilardi, N., Xie, M. H., De Sauvage, F. J., and Gurney, A. L. (2003). Interleukin-23 promotes a distinct CD4 T cell activation state characterized by the production of interleukin-17. J. Biol. Chem. 278, 1910–1914. doi: 10.1074/jbc.M207577200

PubMed Abstract | CrossRef Full Text | Google Scholar

Aguado, E., Richelme, S., Nunez-Cruz, S., Miazek, A., Mura, A. M., Richelme, M., et al. (2002). Induction of T helper type 2 immunity by a point mutation in the LAT adaptor. Science 296, 2036–2040. doi: 10.1126/science.1069057

PubMed Abstract | CrossRef Full Text | Google Scholar

Albert, R., and Thakar, J. (2014). Boolean modeling: a logic-based dynamic approach for understanding signaling and regulatory networks and for making useful predictions. Wiley Interdiscip. Rev. Syst. Biol. Med. 6, 353–369. doi: 10.1002/wsbm.1273

PubMed Abstract | CrossRef Full Text | Google Scholar

Altin, J. A., Tian, L., Liston, A., Bertram, E. M., Goodnow, C. C., and Cook, M. C. (2011). Decreased T cell receptor signaling through CARD11 differentially compromises forkhead box protein 3-positive regulatory versus TH2 effector cells to cause allergy. J. Allergy Clin. Immunol. 127, 1277–1285. doi: 10.1016/j.jaci.2010.12.1081

PubMed Abstract | CrossRef Full Text | Google Scholar

Balázsi, G., Van Oudenaarden, A., and Collins, J. J. (2011). Cellular decision making and biological noise: from microbes to mammals. Cell 144, 910–925. doi: 10.1016/j.cell.2011.01.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Barberis, M., Helikar, T., and Verbruggen, P. (2018). Simulation of stimulation: cytokine dosage and cell cycle crosstalk driving timing-dependent T cell differentiation. Front. Physiol. 9:879. doi: 10.3389/fphys.2018.00879

CrossRef Full Text | Google Scholar

Barberis, M., Todd, R. G., and van der Zee, L. (2017). Advances and challenges in logical modeling of cell cycle regulation: perspective for multi-scale, integrative yeast cell models. FEMS Yeast Res. 17:fow103. doi: 10.1093/femsyr/fow103

PubMed Abstract | CrossRef Full Text | Google Scholar

Barberis, M., and Verbruggen, P. (2017). Quantitative systems biology to decipher design principles of a dynamic cell cycle network: the “Maximum Allowable mammalian Trade–Off–Weight” (MAmTOW). NPJ Syst. Biol. Appl. 3:26. doi: 10.1038/s41540-017-0028-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Becskei, A., and Grusby, M. J. (2007). Contribution of IL-12R mediated feedback loop to Th1 cell differentiation. FEBS Lett. 581, 5199–5206. doi: 10.1016/j.febslet.2007.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Bock, C. N., Babu, S., Breloer, M., Rajamanickam, A., Boothra, Y., Brunn, M.-L., et al. (2017). Th2/1 hybrid cells occurring in murine and human strongyloidiasis share effector functions of Th1 cells. Front. Cell. Infect. Microbiol. 7:261. doi: 10.3389/fcimb.2017.00261

PubMed Abstract | CrossRef Full Text | Google Scholar

Breitfeld, D., Ohl, L., Kremmer, E., Ellwart, J., Sallusto, F., Lipp, M., et al. (2000). Follicular B helper T cells express CXC chemokine receptor 5, localize to B cell follicles, and support immunoglobulin production. J. Exp. Med. 192, 1545–1552. doi: 10.1084/jem.192.11.1545

PubMed Abstract | CrossRef Full Text | Google Scholar

Carbo, A., Hontecillas, R., Kronsteiner, B., Viladomiu, M., Pedragosa, M., Lu, P., et al. (2013). Systems modeling of molecular mechanisms controlling cytokine-driven CD4 + T cell differentiation and phenotype plasticity. PLoS Comput. Biol. 9:e1003027. doi: 10.1371/journal.pcbi.1003027

PubMed Abstract | CrossRef Full Text | Google Scholar

Carbo, A., Olivares-Villagómez, D., Hontecillas, R., Bassaganya-Riera, J., Chaturvedi, R., Blanca Piazuelo, M., et al. (2014). Systems modeling of the role of interleukin-21 in the maintenance of effector CD4 + T cell responses during chronic Helicobacter pylori infection. mBio 5, 1–11. doi: 10.1128/mBio.01243-14

PubMed Abstract | CrossRef Full Text | Google Scholar

Chaouiya, C., Bérenguier, D., Keating, S. M., Naldi, A., van Iersel, M. P., Rodriguez, N., et al. (2013). SBML qualitative models: a model representation format and infrastructure to foster interactions between qualitative modelling formalisms and tools. BMC Syst. Biol. 7:135. doi: 10.1186/1752-0509-7-135

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, L., and Flies, D. B. (2013). Molecular mechanisms of T cell co-stimulation and co-inhibition. Nat. Rev. Immunol. 13, 227–242. doi: 10.1038/nri3405

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W., Jin, W., Hardegen, N., Lei, K., Li, L., Marinos, N., et al. (2003). Conversion of peripheral CD4⁺CD25^-naive T cells to CD4⁺CD25⁺regulatory T cells by TGF-β induction of transcription factor Foxp3. J. Exp. Med. 198, 1875–1886. doi: 10.1084/jem.20030152

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Z., Ding, J., Pang, N., Du, R., and Meng, W. (2013). The Th17/Treg balance and the expression of related cytokines in Uygur cervical cancer patients. Diagn. Pathol. 8:61. doi: 10.1186/1746-1596-8-61

PubMed Abstract | CrossRef Full Text | Google Scholar

Dardalhon, V., Awasthi, A., Kwon, H., Galileos, G., Gao, W., Sobel, R. A., et al. (2008). IL-4 inhibits TGF-beta-induced Foxp3⁺ T cells and, together with TGF-beta, generates IL-9⁺IL-10⁺ Foxp3(-) effector T cells. Nat. Immunol. 9, 1347–1355. doi: 10.1038/ni.1677

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, M.-C., Han, W., Jin, P.-W., Wei, Y.-P., Wei, Q., Zhang, L.-M., et al. (2015). Disturbed Th17/Treg balance in patients with non-small cell lung cancer. Inflammation 38, 2156–2165. doi: 10.1007/s10753-015-0198-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisenstein, E. M., and Williams, C. B. (2009). The Treg/Th17 cell balance: a new paradigm for autoimmunity. Pediatr. Res. 65, 26R–31R. doi: 10.1203/PDR.0b013e31819e76c7

PubMed Abstract | CrossRef Full Text | Google Scholar

Eizenberg-Magar, I., Rimer, J., Zaretsky, I., Lara-Astiaso, D., Reich-Zeliger, S., and Friedman, N. (2017). Diverse continuum of CD4⁺ T cell states is determined by hierarchical additive integration of cytokine signals. Proc. Natl. Acad. Sci. U.S.A. 114, E6447–E6456. doi: 10.1073/pnas.1615590114

PubMed Abstract | CrossRef Full Text | Google Scholar

Evans, C. M., and Jenner, R. G. (2013). Transcription factor interplay in t helper cell differentiation. Brief. Funct. Genomics 12, 499–511. doi: 10.1093/bfgp/elt025

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, D., and Zhu, J. (2017). Dynamic balance between master transcription factors determines the fates and functions of CD4 T cell and innate lymphoid cell subsets. J. Exp. Med. 214, 1861–1876. doi: 10.1084/jem.20170494

PubMed Abstract | CrossRef Full Text | Google Scholar

Gosmann, C., Frazer, I. H., Mattarollo, S. R., and Blumenthal, A. (2014). IL-18, but not IL-12, induces production of IFN-γ in the immunosuppressive environment of HPV16 E7 transgenic hyperplastic skin. J. Invest. Dermatol. 134, 2562–2569. doi: 10.1038/jid.2014.201

PubMed Abstract | CrossRef Full Text | Google Scholar

Groux, H., O’Garra, A., Bigler, M., Rouleau, M., Antonenko, S., de Vries, J. E., et al. (1997). A CD4⁺ T cell subset inhibits antigen-specific T cell responses and prevents colitis. Nature 389, 737–742. doi: 10.1038/39614

PubMed Abstract | CrossRef Full Text | Google Scholar

Harbour, S. N., Maynard, C. L., Zindl, C. L., Schoeb, T. R., and Weaver, C. T. (2015). Th17 cells give rise to Th1 cells that are required for the pathogenesis of colitis. Proc. Natl. Acad. Sci. U.S.A. 112, 7061–7066. doi: 10.1073/pnas.1415675112

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrington, L. E., Hatton, R. D., Mangan, P. R., Turner, H., Murphy, T. L., Murphy, K. M., et al. (2005). Interleukin 17–producing CD4⁺ effector T cells develop via a lineage distinct from the T helper type 1 and 2 lineages. Nat. Immunol. 6, 1123–1132. doi: 10.1038/ni1254

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrison, O. J., Srinivasan, N., Pott, J., Schiering, C., Krausgruber, T., Ilott, N. E., et al. (2015). Epithelial-derived IL-18 regulates Th17 cell differentiation and Foxp3⁺ Treg cell function in the intestine. Mucosal Immunol. 8, 1226–1236. doi: 10.1038/mi.2015.13

PubMed Abstract | CrossRef Full Text | Google Scholar

Hegazy, A. N., Peine, M., Helmstetter, C., Panse, I., Fröhlich, A., Bergthaler, A., et al. (2010). Interferons direct Th2 cell reprogramming to generate a stable GATA-3⁺ T-bet + Cell Subset with Combined Th2 and Th1 Cell Functions. Immunity 32, 116–128. doi: 10.1016/j.immuni.2009.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Helikar, T., Konvalina, J., Heidel, J., and Rogers, J. A. (2008). Emergent decision-making in biological signal transduction networks. Proc. Natl. Acad. Sci. U.S.A. 105:705088105. doi: 10.1073/pnas.0705088105

PubMed Abstract | CrossRef Full Text | Google Scholar

Helikar, T., Kowal, B., Madrahimov, A., Shrestha, M., Pedersen, J., Limbu, K., et al. (2012a). Bio-logic builder: a Non-Technical tool for building dynamical, qualitative models. PLoS One 7:e46417. doi: 10.1371/journal.pone.0046417

PubMed Abstract | CrossRef Full Text | Google Scholar

Helikar, T., Kowal, B., McClenathan, S., Bruckner, M., Rowley, T., Madrahimov, A., et al. (2012b). The cell collective: toward an open and collaborative approach to systems biology. BMC Syst. Biol. 6:96. doi: 10.1186/1752-0509-6-96

PubMed Abstract | CrossRef Full Text | Google Scholar

Helikar, T., Kowal, B., and Rogers, J. A. (2013). A cell simulator platform: the cell collective. Clin. Pharmacol. Ther. 93, 393–395. doi: 10.1038/clpt.2013.41

PubMed Abstract | CrossRef Full Text | Google Scholar

Helikar, T., and Rogers, J. A. (2009). ChemChains: a platform for simulation and analysis of biochemical networks aimed to laboratory scientists. BMC Syst. Biol. 3:58. doi: 10.1186/1752-0509-3-58

PubMed Abstract | CrossRef Full Text | Google Scholar

Ichiyama, K., Yoshida, H., Wakabayashi, Y., Chinen, T., Saeki, K., Nakaya, M., et al. (2008). Foxp3 inhibits RORgammat-mediated IL-17A mRNA transcription through direct interaction with RORgammat. J. Biol. Chem. 283, 17003–17008. doi: 10.1074/jbc.M801286200

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaiko, G. E., Horvat, J. C., Beagley, K. W., and Hansbro, P. M. (2008). Immunological decision-making: how does the immune system decide to mount a helper T cell response? Immunology 123, 326–338. doi: 10.1111/j.1365-2567.2007.02719.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaplan, M. H. (2013). Th9 cells: differentiation and disease. Immunol. Rev. 252, 104–115. doi: 10.1111/imr.12028

PubMed Abstract | CrossRef Full Text | Google Scholar

Kimura, A., and Kishimoto, T. (2010). IL-6: regulator of Treg/Th17 balance. Eur. J. Immunol. 40, 1830–1835. doi: 10.1002/eji.201040391

PubMed Abstract | CrossRef Full Text | Google Scholar

Koch, M. A., Tucker-Heard, G., Perdue, N. R., Killebrew, J. R., Urdahl, K. B., and Campbell, D. J. (2009). The transcription factor T-bet controls regulatory T cell homeostasis and function during type 1 inflammation. Nat. Immunol. 10, 595–602. doi: 10.1038/ni.1731

PubMed Abstract | CrossRef Full Text | Google Scholar

Kullberg, M. C., Jankovic, D., Feng, C. G., Hue, S., Gorelick, P. L., McKenzie, B. S., et al. (2006). IL-23 plays a key role in Helicobacter hepaticus-induced T cell-dependent colitis. J. Exp. Med. 203, 2485–2494. doi: 10.1084/jem.20061082

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Novère, N. (2015). Quantitative and logic modelling of molecular and gene networks. Nat. Rev. Genet. 16, 146–158. doi: 10.1038/nrg3885

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, Y. K., Mukasa, R., Hatton, R. D., and Weaver, C. T. (2009a). Developmental plasticity of Th17 and Treg cells. Curr. Opin. Immunol. 21, 274–280. doi: 10.1016/j.coi.2009.05.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, Y. K., Turner, H., Maynard, C. L., Oliver, J. R., Chen, D., Elson, C. O., et al. (2009b). Late Developmental Plasticity in the T Helper 17 Lineage. Immunity 30, 92–107. doi: 10.1016/j.immuni.2008.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Linke, C., Chasapi, A., González-Novo, A., Al Sawad, I., Tognetti, S., Klipp, E., et al. (2017). A Clb/Cdk1-mediated regulation of Fkh2 synchronizes CLB expression in the budding yeast cell cycle. NPJ Syst. Biol. Appl. 3:7. doi: 10.1038/s41540-017-0008-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Luchting, B., Rachinger-Adam, B., Zeitler, J., Egenberger, L., Möhnle, P., Kreth, S., et al. (2014). Disrupted TH17/Treg balance in patients with chronic low back pain. PLoS ONE 9:e104883. doi: 10.1371/journal.pone.0104883

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, C. S., Deenick, E. K., Batten, M., and Tangye, S. G. (2012). The origins, function, and regulation of T follicular helper cells. J. Exp. Med. 209, 1241–1253. doi: 10.1084/jem.20120994

PubMed Abstract | CrossRef Full Text | Google Scholar

Mangan, P. R., Harrington, L. E., O’Quinn, D. B., Helms, W. S., Bullard, D. C., Elson, C. O., et al. (2006). Transforming growth factor-β induces development of the TH17 lineage. Nature 441, 231–234. doi: 10.1038/nature04754

PubMed Abstract | CrossRef Full Text | Google Scholar

Marino, S., Hogue, I. B., Ray, C. J., and Kirschner, D. E. (2008). A methodology for performing global uncertainty and sensitivity analysis in systems biology. J. Theor. Biol. 254, 178–196. doi: 10.1016/j.jtbi.2008.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Martinez-Sanchez, M. E., Huerta, L., Alvarez-Buylla, E. R., and Villarreal Lujan, C. (2018). Role of cytokine combinations on CD4+ T cell differentiation, partial polarization, and plasticity: continuous network modeling approach. Front. Physiol. 9:877. doi: 10.3389/fphys.2018.00877

CrossRef Full Text | Google Scholar

Martinez-Sanchez, M. E., Mendoza, L., Villarreal, C., and Alvarez-Buylla, E. R. (2015). A Minimal Regulatory Network of Extrinsic and Intrinsic Factors Recovers Observed Patterns of CD4 + T Cell Differentiation and Plasticity. PLoS Comput. Biol. 11:e1004324. doi: 10.1371/journal.pcbi.1004324

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrison, P. J., Bending, D., Fouser, L. A., Wright, J. F., Stockinger, B., Cooke, A., et al. (2013). Th17-cell plasticity in Helicobacter hepaticus-induced intestinal inflammation. Mucosal Immunol. 6, 1143–1156. doi: 10.1038/mi.2013.11

PubMed Abstract | CrossRef Full Text | Google Scholar

Mosmann, T. R., Cherwinski, H., Bond, M. W., Giedlin, M. A., and Coffman, R. L. (1986). Two types of murine helper T cell clone. I. Definition according to profiles of lymphokine activities and secreted proteins. J. Immunol. 136, 2348–2357. doi: 10.1111/j.1442-9071.2011.02672.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Munk, R. B., Sugiyama, K., Ghosh, P., Sasaki, C. Y., Rezanka, L., Banerjee, K., et al. (2011). Antigen-independent IFN-γ production by human naïve CD4⁺ T cells activated by IL-12 plus IL-18. PLoS One 6:e18553. doi: 10.1371/journal.pone.0018553

PubMed Abstract | CrossRef Full Text

Murai, H., Qi, H., Choudhury, B., Wild, J., Dharajiya, N., Vaidya, S., et al. (2012). Alternaria-induced release of IL-18 from damaged airway epithelial cells: an NF-??B dependent mechanism of Th2 differentiation? PLoS One 7:e30280. doi: 10.1371/journal.pone.0030280

PubMed Abstract | CrossRef Full Text | Google Scholar

Murphy, K. M., and Reiner, S. L. (2002). The lineage decisions of helper T cells. Nat. Rev. Immunol. 2, 933–944. doi: 10.1038/nri954

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakanishi, K., Yoshimoto, T., Tsutsui, H., and Okamura, H. (2001). Interleukin-18 regulates both Th1 and Th2 responses. Annu. Rev. Immunol. 19, 423–474. doi: 10.1146/annurev.immunol.19.1.423

CrossRef Full Text | Google Scholar

Naldi, A., Carneiro, J., Chaouiya, C., and Thieffry, D. (2010). Diversity and plasticity of Th cell types predicted from regulatory network modelling. PLoS Comput. Biol. 6:e1000912. doi: 10.1371/journal.pcbi.1000912

PubMed Abstract | CrossRef Full Text | Google Scholar

Naldi, A., Monteiro, P. T., Müssel, C., Kestler, H. A., Thieffry, D., Xenarios, I., et al. (2015). Cooperative development of logical modelling standards and tools with CoLoMoTo. Bioinformatics 31, 1154–1159. doi: 10.1093/bioinformatics/btv013

PubMed Abstract | CrossRef Full Text | Google Scholar

Nindl, V., Maier, R., Ratering, D., De Giuli, R., Züst, R., Thiel, V., et al. (2012). Cooperation of Th1 and Th17 cells determines transition from autoimmune myocarditis to dilated cardiomyopathy. Eur. J. Immunol. 42, 2311–2321. doi: 10.1002/eji.201142209

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohnmacht, C., Park, J., Cording, S., Wing, J. B., Atarashi, K., Obata, Y., et al. (2015). The microbiota regulates type 2 immunity through RORγt + T cells. Science 349, 1–9. doi: 10.1126/science.aac4263

PubMed Abstract | CrossRef Full Text | Google Scholar

Oldenhove, G., Bouladoux, N., Wohlfert, E. A., Hall, J. A., Chou, D., Dos Santos, L., et al. (2009). Decrease of Foxp3 + Treg cell number and acquisition of effector cell phenotype during lethal infection. Immunity 31, 772–786. doi: 10.1016/j.immuni.2009.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Omenetti, S., and Pizarro, T. T. (2015). The Treg/Th17 axis: a dynamic balance regulated by the gut microbiome. Front. Immunol. 6:639. doi: 10.3389/fimmu.2015.00639

PubMed Abstract | CrossRef Full Text | Google Scholar

Palau-Ortin, D., Formosa-Jordan, P., Sancho, J. M., and Ibañes, M. (2015). Pattern selection by dynamical biochemical signals. Biophys. J. 108, 1555–1565. doi: 10.1016/j.bpj.2014.12.058

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, H., Li, Z., Yang, X. O., Chang, S. H., Nurieva, R., Wang, Y.-H., et al. (2005). A distinct lineage of CD4 T cells regulates tissue inflammation by producing interleukin 17. Nat. Immunol. 6, 1133–1141. doi: 10.1038/ni1261

PubMed Abstract | CrossRef Full Text | Google Scholar

Paul, W. E. (2010). What determines Th2 differentiation, in vitro and in vivo? Immunol. Cell Biol. 88, 236–239. doi: 10.1038/icb.2010.2

PubMed Abstract | CrossRef Full Text | Google Scholar

Peine, M., Rausch, S., Helmstetter, C., Fröhlich, A., Hegazy, A. N., Kühl, A. A., et al. (2013). Stable T-bet + GATA-3 + Th1/Th2 hybrid cells arise in vivo, can develop directly from naive precursors, and limit immunopathologic inflammation. PLoS Biol. 11:e1001633. doi: 10.1371/journal.pbio.1001633

PubMed Abstract | CrossRef Full Text | Google Scholar

Podojil, J. R., and Miller, S. D. (2009). Molecular mechanisms of T cell receptor and costimulatory molecule ligation/blockade in autoimmune disease therapy. Immunol. Rev. 229, 337–355. doi: 10.1111/j.1600-065X.2009.00773.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Pujol, G., Iooss, B., Alexandre Janon with contributions from Khalid Boumhaout, Da Veiga, S. A., Delage, T., Fruth, J., et al. (2017). Sensitivity: Global Sensitivity Analysis of Model Outputs. Available at: https://cran.r-project.org/package=sensitivity

Google Scholar

R Development Core Team (2011). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing, doi: 10.1007/978-3-540-74686-7

CrossRef Full Text | Google Scholar

Reiner, S. L. (2007). Development in motion: helper T cells at work. Cell 129, 33–36. doi: 10.1016/j.cell.2007.03.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Ribeiro, A. S., and Kauffman, S. A. (2007). Noisy attractors and ergodic sets in models of gene regulatory networks. J. Theor. Biol. 247, 743–755. doi: 10.1016/j.jtbi.2007.04.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Romagnani, S. (2000). T cell subsets (Th1 versus Th2). Ann. Allergy Asthma Immunol. 85, 9–18; quiz 18, 21. doi: 10.1016/S1081-1206(10)62426-X

CrossRef Full Text | Google Scholar

Rowell, E., and Wilson, C. B. (2009). Programming perpetual T helper cell plasticity. Immunity 30, 7–9. doi: 10.1016/j.immuni.2008.12.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Samaga, R., and Klamt, S. (2013). Modeling approaches for qualitative and semi-quantitative analysis of cellular signaling networks. Cell Commun. Signal. 11:43. doi: 10.1186/1478-811X-11-43

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaerli, P., Willimann, K., Lang, A. B., Lipp, M., Loetscher, P., and Moser, B. (2000). CXC chemokine receptor 5 expression defines follicular homing T cells with B cell helper function. J. Exp. Med. 192, 1553–1562. doi: 10.1084/jem.192.11.1553

CrossRef Full Text | Google Scholar

Schmitt, E. G., and Williams, C. B. (2013). Generation and function of induced regulatory T cells. Front. Immunol. 4:152. doi: 10.3389/fimmu.2013.00152

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, G., Cox, C. A., Vistica, B. P., Tan, C., Wawrousek, E. F., and Gery, I. (2008). Phenotype switching by inflammation-inducing polarized Th17 cells, but not by Th1 cells. J. Immunol. 181, 7205–7213. doi: 10.4049/jimmunol.181.10.7205

PubMed Abstract | CrossRef Full Text | Google Scholar

Soroosh, P., and Doherty, T. A. (2009). Th9 and allergic disease. Immunology 127, 450–458. doi: 10.1111/j.1365-2567.2009.03114.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tarjan, R. (1972). Depth-first search and linear graph algorithms. SIAM J. Comput. 1, 146–160. doi: 10.1137/0201010

PubMed Abstract | CrossRef Full Text | Google Scholar

Tartar, D. M., VanMorlan, A. M., Wan, X., Guloglu, F. B., Jain, R., Haymaker, C. L., et al. (2010). FoxP3⁺ ROR t + T helper intermediates display suppressive function against autoimmune diabetes. J. Immunol. 184, 3377–3385. doi: 10.4049/jimmunol.0903324

CrossRef Full Text | Google Scholar

Tesmer, L. A., Lundy, S. K., Sarkar, S., and Fox, D. A. (2008). Th17 cells in human disease. Immunol. Rev. 223, 87–113. doi: 10.1111/j.1600-065X.2008.00628.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Tijms, H. C. (2003). A First Course in Stochastic Models. New York, NY: John Wiley & Sons, Ltd, doi: 10.1002/047001363X

CrossRef Full Text | Google Scholar

Todd, R. G., and Helikar, T. (2012). Ergodic sets as cell phenotype of budding yeast cell cycle. PLoS One 7:e45780. doi: 10.1371/journal.pone.0045780

PubMed Abstract | CrossRef Full Text | Google Scholar

Tominaga, K., Yoshimoto, T., Torigoe, K., Kurimoto, M., Matsui, K., Hada, T., et al. (2000). IL-12 synergizes with IL-18 or IL-1beta for IFN-gamma production from human T cells. Int. Immunol. 12, 151–160. doi: 10.1093/intimm/12.2.151

CrossRef Full Text | Google Scholar

van Panhuys, N., Klauschen, F., and Germain, R. N. (2014). T cell-receptor-dependent signal intensity dominantly controls CD4 + T cell polarization in vivo. Immunity 41, 63–74. doi: 10.1016/j.immuni.2014.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Veldhoen, M., Uyttenhove, C., van Snick, J., Helmby, H., Westendorf, A., Buer, J., et al. (2008). Transforming growth factor-beta “reprograms” the differentiation of T helper 2 cells and promotes an interleukin 9-producing subset. Nat. Immunol. 9, 1341–1346. doi: 10.1038/ni.1659

PubMed Abstract | CrossRef Full Text | Google Scholar

Voo, K. S., Wang, Y.-H., Santori, F. R., Boggiano, C., Wang, Y.-H., Arima, K., et al. (2009). Identification of IL-17-producing FOXP3 + regulatory T cells in humans. Proc. Natl. Acad. Sci. U.S.A. 106, 4793–4798. doi: 10.1073/pnas.0900408106

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, L., Kitani, A., Fuss, I., and Strober, W. (2007). Cutting edge: regulatory T cells induce CD4⁺ CD25-Foxp3- T cells or are self-induced to become Th17 cells in the absence of exogenous TGF-beta. J. Immunol. 178, 6725–6729. doi: 10.4049/jimmunol.178.11.6725

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamane, H., and Paul, W. E. (2012). Cytokines of the γc family control CD4 + T cell differentiation and function. Nat. Immunol. 13, 1037–1044. doi: 10.1038/ni.2431

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamane, H., Zhu, J., and Paul, W. E. (2005). Independent roles for IL-2 and GATA-3 in stimulating naive CD4 + T cells to generate a Th2-inducing cytokine environment. J. Exp. Med. 202, 793–804. doi: 10.1084/jem.20051304

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X. O., Nurieva, R., Martinez, G. J., Kang, H. S., Chung, Y., Pappu, B. P., et al. (2008). Molecular antagonism and plasticity of regulatory and inflammatory T cell programs. Immunity 29, 44–56. doi: 10.1016/j.immuni.2008.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshimoto, T., Takeda, K., Tanaka, T., Ohkusu, K., Kashiwamura, S., Okamura, H., et al. (1998). IL-12 up-regulates IL-18 receptor expression on T cells, Th1 cells, and B cells: synergism with IL-18 for IFN-gamma production. J. Immunol. 161, 3400–3407.

PubMed Abstract | Google Scholar

Yu, F., Sharma, S., Edwards, J., Feigenbaum, L., and Zhu, J. (2015). Dynamic expression of transcription factors T-bet and GATA-3 by regulatory T cells maintains immunotolerance. Nat. Immunol. 16, 197–206. doi: 10.1038/ni.3053

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, L., Lopes, J. E., Chong, M. M. W., Ivanov, I. I., Min, R., Victora, G. D., et al. (2008). TGF-β-induced Foxp3 inhibits TH17 cell differentiation by antagonizing RORγt function. Nature 453, 236–240. doi: 10.1038/nature06878

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, J., and Paul, W. E. (2008). CD4 T cells: fates, functions, and faults. Blood 112, 1557–1569. doi: 10.1182/blood-2008-05-078154

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, J., Yamane, H., and Paul, W. E. (2010). Differentiation of effector CD4 T cell populations. Annu. Rev. Immunol. 28, 445–489. doi: 10.1146/annurev-immunol-030409-101212

PubMed Abstract | CrossRef Full Text | Google Scholar