Although the evidence is still limited, a growing body of research suggests music may have beneficial effects for diseases such as Parkinson’s.
Image credit: Shutterstock/agsandrew.
New Research In
Physical Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
Edited by Federica Sallusto, Institute for Research in Biomedicine, Università della Svizzera Italiana, Bellinzona, Switzerland, and accepted by Editorial Board Member Tadatsugu Taniguchi June 2, 2017 (received for review September 28, 2016)
Significance
Understanding the logic by which cells respond to complex signal combinations is challenging. We used CD4+ T cells as a model system to study signal integration by systematically mapping their differentiation in response to a large number of cytokine combinations. We find that, in response to varied cytokine mixtures, cells coexpress lineage-specifying proteins at diverse levels, such that the cell population spans a continuum of intermediate states between canonical cell phenotypes. Mathematical modeling explains these results using hierarchical summation of cytokine inputs and correctly predicts population response to new input conditions. These findings suggest that complex cellular responses can be effectively described using relatively simple hierarchical summation rules, providing a framework for prediction of cellular responses to signal combinations.
Abstract
During cell differentiation, progenitor cells integrate signals from their environment that guide their development into specialized phenotypes. The ways by which cells respond to complex signal combinations remain difficult to analyze and model. To gain additional insight into signal integration, we systematically mapped the response of CD4+ T cells to a large number of input cytokine combinations that drive their differentiation. We find that, in response to varied input combinations, cells differentiate into a continuum of cell fates as opposed to a limited number of discrete phenotypes. Input cytokines hierarchically influence the cell population, with TGFβ being most dominant followed by IL-6 and IL-4. Mathematical modeling explains these results using additive signal integration within hierarchical groups of input cytokine combinations and correctly predicts cell population response to new input conditions. These findings suggest that complex cellular responses can be effectively described using a segmented linear approach, providing a framework for prediction of cellular responses to new cytokine combinations and doses, with implications to fine-tuned immunotherapies.
Cell differentiation is controlled by complex gene regulatory networks that determine the expression of a large number of genes in response to external stimuli. CD4+ T cells, central regulators of immune responses, can differentiate into a number of phenotypes (or cell states), each defined by the expression of a panel of genes, such as lineage-specifying transcription factors (TFs) and cytokines (1). Differentiation outcome depends on a number of factors, including strength of antigen stimulation (2, 3), interactions with antigen presenting cells, and the signaling environment defined by secreted cytokines (1). By applying specific cytokine signals, it is possible to direct antigen-activated T cells toward differentiation into a number of specific phenotypes. Although the molecular interactions and regulatory networks that govern the differentiation process are being revealed at increasing resolution (4, 5), understanding the logic of operation of these complex networks and developing predictive mathematical models for cellular responses remain a challenge.
Mathematical models that describe complex cellular processes are difficult to construct because of a lack of complete information about the underlying networks of molecular interactions and in particular, missing quantitative data regarding the parameters of biochemical interactions within these networks. Alternatively, a “black box” approach can be used, constructing a model of the system by analyzing its response to different input conditions (Fig. 1A). This technique enables study of complex systems in a simplified manner, because its implementation merely requires measurable inputs and outputs, without quantifying the response of each component inside the box. Approaches based on this concept are highly useful for describing engineered systems and were also applied for studying some biological systems, such as signal transduction in the osmoadaptation pathway of yeast (6) and agent-based modeling of multicellular interactions (7). Here, we apply the black box approach to construct a predictive model of cell differentiation. We use CD4+ T-cell differentiation as a model system and construct a mathematical model that describes cell differentiation based on measurements of cell phenotype in response to a large number of input cytokines and their combinations.
Fig. 1.
Systematic characterization of input–output relations in CD4+ T-cell differentiation. (A) Schematic representation of the black box concept. Cells are cultured with combinations of input signals (I1, I2,…, IN), and their phenotype (cell state) is determined by measuring expression levels of a set of lineage characteristic proteins (output: O1, O2,…, OM) without determining levels of all nodes in the regulatory network that controls cellular responses. (B) Experimental design. Naïve CD4+ T cells were provided with TCR stimulation (anti-CD3) and costimulation (anti-CD28) and cultured with all possible 64 combinations of six input cytokines. Each cytokine was either used at a fixed concentration or absent from the cell culture medium. (C) After 7 d of culture, cells were restimulated, and levels of 10 proteins that are characteristic of cell differentiation state were measured using flow cytometry. Examples of protein-level distributions are shown for four input combinations. Input cytokine combinations are marked on the left side using the one-symbol labels as in B (β, TGFβ; γ, IFNγ; t, IL-12; 2, IL-2; 4, IL-4; 6, IL-6).
CD4+ T cells are a major hub of immune regulation. Depending on signals sensed from their microenvironment during initial antigen priming, naïve CD4+ T cells can adopt various effector and regulatory phenotypes, including Th1, Th2, Th17, and induced Treg (iTreg) cells (1, 8, 9). These distinct phenotypes differ from one another in their genomic and epigenomic profile and their functional characteristics. Th1 cells require IL-12 to induce their development and are characterized by the expression of IFNγ and the lineage-specifying TF T-bet (10⇓–12). IL-4 presence, further supported by IL-2 (13, 14), drives the differentiation of Th2 cells, which are characterized by the induction of the key regulator GATA3 and the production of the cytokines IL-4, IL-5, and IL-13 (15⇓–17). Th17 cell differentiation is induced in the presence of TGFβ together with the proinflammatory cytokine IL-6 and is characterized by the production of IL-17A, IL-17F, IL-21, and IL-22, and the expression of the TF RORγt (18⇓⇓⇓–22). Treg cells are characterized by expression of the TF Foxp3, which can be induced in peripheral naïve CD4+ T cells by a combination of TGFβ and IL-2. Foxp3 drives the expression of the cytokines TGFβ, IL-10, and IL-35, which play major roles in Treg suppressive activity (23⇓⇓⇓⇓–28). Functionally, these diversified CD4+ T-cell subsets coordinate different sets of immune responses required to protect the body against various infections; unregulated or insufficient response may result in inflammatory or autoimmune pathogenesis.
Common perception held that each CD4+ T-cell lineage requires different inducing conditions and possesses a distinct cytokine profile. However, this notion is challenged by evidence of population heterogeneity (29⇓⇓–32) and Th cell plasticity (33⇓⇓⇓–37), showing that, under various conditions, differentiated cells can adopt hybrid phenotypes, coexpressing signature TFs and cytokines of different lineages. We and others have recently shown that naïve CD4+ T cells can differentiate into a mixed Th1–Th2 phenotype, coexpressing lineage-specific TFs and cytokines of both cell types in varying levels depending on the relative levels of input cytokine signals (38⇓–40). Taken together, these observations suggest a more complex picture of CD4+ T-cell fate determination with a substantial degree of phenotype plasticity and flexibility in response to different input stimuli.
Although CD4+ T-cell plasticity is gaining appreciation as a central characteristic of their function (41), the logic by which different input signals are combined to drive complex cell states remains obscure. Thus, we set out to investigate how CD4+ T cells respond to a number of simultaneous signals and combine them to define a specific differentiation program. To gain a systematic view of signal integration, we explored the differentiation of CD4+ T cells in an unbiased way by exposing cells to a large number of cytokine signal combinations and evaluating the response of the cell population. Specifically, we used combinations of six input cytokines that drive the differentiation of the four “classical” CD4+ T-cell lineages described above (Th1, Th2, Th17, and iTreg). After culture, we measured the level of expression of 10 lineage-specifying TFs and cytokines (output; see below). We use the levels of these 10 proteins to define cells’ state in a 10D “differentiation space.” Our study revealed a large number of possible cell states, in which cells express diverse combinations of cytokines and TFs at varied intermediate expression levels.
Although the response of individual cells is heterogeneous, with considerable cell to cell variability among cells cultured under the same conditions, the response of the cell population is more predictable. The population-averaged cell states measured under the different input combinations did not cluster into a small number of distinct phenotypes but were scattered in a large region of the differentiation space, forming a continuum that spans the range between the four classical phenotypes. Organization of these mixed cell states in the differentiation space reflected a hierarchy in the influence of the different cytokines on cell differentiation, with TGFβ being the most potent input. Based on these observations, we constructed a mathematical model describing the differentiation process, in which cell populations respond by additively integrating signals from the combinatorial cytokine milieu in a hierarchical manner. The model can predict the response of the cell population to new input conditions, which was validated experimentally.
Results
Mapping CD4+ T-Cell Differentiation Under a Large Number of Mixed Input Conditions Reveals Heterogeneous Cellular Responses.
We mapped the differentiation decision space of CD4+ T cells by examining the phenotypes of cells cultured with a large number of combinations of cytokines. Naïve CD4+ T cells were isolated from spleens of C57BL/6 mice and cultured with stimulation of their T-cell receptor (TCR) and costimulation using plate-bound αCD3 and αCD28 antibodies. The cells were supplemented with combinations of six cytokines (IFNγ, IL-12, TGFβ, IL-6, IL-4, and IL-2), which are known inducers of four main Th subsets (Th1, Th2, Th17, and iTreg) (1). In these experiments, each of six input cytokines was used at a fixed concentration, which we found in calibration experiments to induce the differentiation of cells toward these four Th lineages (SI Appendix, Fig. S1). We used all binary combinations of these six cytokines, with each cytokine either present or absent: a total of 26 = 64 conditions (Fig. 1B). It should be noted that all of the cytokines used as inputs, other than IL-12, can also be secreted by CD4+ T cells following activation and therefore, can be present in the culture medium even under conditions in which they were not externally added.
After 7 d of culture, cells were restimulated, and their cytokine and TF expression levels were measured using intracellular antibody staining followed by flow cytometry. Substantial cell proliferation was observed under all mixed conditions, with relatively small differences in proliferation capacity between different input mixtures (SI Appendix, Fig. S2). We obtained 10 different fluorescence readings characterizing cell state under each of 64 input conditions: six output cytokines (IL-4, IL-5, IL-13, IFNγ, IL-10, and IL-17A) and four lineage-specifying TFs (GATA3, T-bet, Foxp3, and RORγt). For analysis of the population-averaged response, we used normalized mean fluorescence intensity (MFI) values of the measured outputs (details are in Experimental Procedures).
Fig. 1C shows the measured levels of 10 proteins for four representative input conditions. As expected, adding IL-4 + IL-2 (Fig. 1C, row 1) to the cell culture resulted in high expression of the Th2-associated cytokines IL-4, IL-5, and IL-13 and the TF GATA3. IL-10 level is also high under this condition, in line with previous observations (42), whereas levels of the other proteins are low or intermediate. Adding IL-12 to the medium resulted in high expression of IFNγ and T-bet (Fig. 1C, row 2), which are associated with Th1, and relatively low levels of the other proteins tested. Adding TGFβ + IL-2 resulted in a bimodal expression of Foxp3, with the majority of cells expressing this TF, as expected (Fig. 1C, row 3). However, the global mapping of responses to 64 cytokine combinations reveals a more complex picture. For example, in the presence of IFNγ + IL-12 + IL-6 + IL-4 + IL-2, the levels of several proteins, associated with different lineages, were elevated (Fig. 1C, row 4). In general, expected patterns of expression are observed for conditions that are known to provide a defined Th cell phenotype, whereas more complex input combinations resulted in diverse outcomes, with proteins expressed at various combinations and in different intermediate levels.
Mean Expression Levels of Lineage-Specific Proteins Span a Continuum of Intermediate Values.
We continue our analysis by examining the distributions of expression levels for each of 10 proteins measured. All proteins displayed a relatively broad distribution in their expression level across the cell population under most conditions. Examples for the distribution of expression levels of GATA3 and IFNγ under a number of input combinations are shown in Fig. 2 A and B, respectively. Levels of the TF GATA3 were high under conditions that favor Th2 differentiation, such as IL-4 + IL-2. GATA3 levels gradually dropped under different input combinations (such as IL-4 + IL-12 > IL-2 + IL-12 > IL-12), reaching very low levels of expression under conditions, such as TGFβ + IL-6 + IFNγ (Fig. 2A). Similar results were also obtained for the other measured TFs and cytokines (Fig. 2 B–D and SI Appendix, Fig. S3).
Fig. 2.
Continuity in expression levels of lineage-specific cytokines and TFs under mixed input conditions. (A and B) Distributions of GATA3 (A) and IFNγ (B) protein expression levels under different cytokine combinations measured using fluorescently labeled antibodies and flow cytometry. Expression levels are heterogeneous within populations of cells cultured under the same condition, and mean expression levels span a continuum of intermediate values. (C and D) Distributions of T-bet (C) and TNFα (D) protein expression levels under different cytokine combinations measured using mass cytometry. Inset in C compares T-bet mean expression levels measured by fluorescence and mass cytometry. Each dot represents normalized mean levels measured under a specific input combination (SI Appendix, Fig. S3). (E and F) Comparison of mean GATA3 (E) and mean IL-10 (F) mRNA and protein expression levels. Protein levels were measured by antibody labeling and flow cytometry (x axis), and mRNA levels were measured by RNA-seq (y axis). Each dot represents normalized mean levels measured under a specific input combination (Data Normalization and Analysis). (G and H) Mean RORγt expression levels under a conventional Th17-inducing condition (TGFβ + IL-6), a conventional iTreg-inducing condition (TGFβ + IL-2), and two additional cytokine combinations. Adding IL-4, IL-2, and IL-12 supports elevated RORγt levels as measured for both protein (G) and mRNA (H). (I and J) IL-17A expression levels under the same conditions as in G and H for protein (I) and mRNA (J). Expression levels in G–J were normalized by the level measured in the TGFβ + IL-2 sample. Input cytokine combinations are marked using the one-symbol labels as specified in the table at the bottom.
For some of the conditions, we also repeated these measurements using mass cytometry (CyTOF). We obtained similar results in both techniques, with lineage-associated TFs and cytokines showing a range of expression levels in response to various mixed inputs (Fig. 2 C and D and SI Appendix, Fig. S3). We further measured levels of gene expression using RNA sequencing (RNA-seq) for cells cultured under a set of various input combinations. In general, measured mRNA levels correlate well with mean protein levels as measured by flow cytometry (Fig. 2 E and F). Intermediate expression levels are also observed for mRNA under various input mixtures, ranging between high and low levels that were observed for known polarizing conditions.
Our approach for mapping T-cell differentiation under a large number of cytokine mixtures allowed us to study cellular responses to conditions that are not typically investigated. We find that lineage characteristic proteins can be expressed at high levels under conditions that were not previously considered standard for driving their expression. For example, levels of RORγt, a key regulator for Th17 differentiation, were high in the presence of TGFβ + IL-6, as was previously shown (19, 20); however, we find that its expression level is also elevated when other cytokines (IL-4, IL-2, and IL-12) were added to the mixture (Fig. 2G). These input conditions also support the expression of IL-17A, a characteristic cytokine secreted by Th17 cells (Fig. 2I). These observations are supported by the elevated mRNA levels of RORγt and IL-17 that we measured under these input combinations (Fig. 2 H and J). Thus, Th1- and Th2-inducing cytokines are supportive for the expression of Th17-associated proteins rather than suppressive as might be expected.
Together, our results show that high average levels of each lineage-specific protein are obtained under conditions known to drive differentiation toward this lineage (Fig. 2). However, under other cytokine conditions, the expression of these proteins is not necessarily low; rather, the distributions continuously shift between low and high levels, and the expression is intermediate under various input combinations (Fig. 2 and SI Appendix, Fig. S3).
CD4+ T Cells Show Complex Patterns of Coexpression of Lineage-Specific Proteins in Response to Combinations of Input Signals.
Next, we examined protein coexpression patterns at the single-cell level. Coexpression was observed for many of the tested proteins under most experimental conditions, including proteins characterizing different lineages. Conventional polarizing conditions resulted in the expression of hallmark cytokines and TFs, which characterize the different polarized Th lineages (Fig. 3 A, columns 1 and 4 and B, columns 1 and 4). However, other cytokine combinations elicited mixed phenotypes. On stimulation in the presence of IL-12 + IL-2 or IL-12 + IL-4, we detected a subpopulation of cells coexpressing IFNγ, a Th1 cytokine, and IL-4, a Th2 cytokine (Fig. 3A, columns 2 and 3), extending previous results (38). In the presence of IFNγ + TGFβ + IL-6 + IL-2 or TGFβ + IL-4 + IL-2, we find cells coexpressing RORγt, a Th17 regulator, together with Foxp3, a regulator of Treg cells (Fig. 3B, columns 2 and 3). In a similar manner, most cytokine input combinations that we tested drove mixed patterns of coexpression of several proteins, resulting in a heterogeneous cell population. To examine whether the mixed states represent a transient phase adopted by the cells during the course of their differentiation toward polarized lineages, we extended the duration of culture to 14 d. We characterized the expression levels of four lineage-associated TFs (T-bet, GATA3, RORγt, and Foxp3) in cells cultured under four representative mixed input conditions as well as cells cultured under four polarizing conditions. As can be expected, polarizing conditions after 7 and 14 d of culture resulted in high expression of the characteristic lineage-specifying TF and low expression of the other three TFs (Fig. 3C, Left). In contrast, more than one lineage-specifying TF is expressed at intermediate to high levels in each of the tested mixed input conditions at both time points (Fig. 3C, Right). The overall patterns of TFs coexpression remain rather stable between 7 and 14 d of culture. Therefore, doubling the time of culture did not result in a convergence of the mixed phenotypes toward a more polarized state; rather, the population still exhibited a mixed phenotype after the prolonged culture.
Fig. 3.
Patterns of protein coexpression and stability of TF expression under cytokine input combinations. (A) Scatter plots showing coexpression patterns of Th1- and Th2-associated cytokines (IFNγ and IL-4, respectively) under different input combinations: column 1, Th1-inducing conditions; column 2, IL-12 + IL-2; column 3, IL-12 + IL-4; column 4, Th2-inducing conditions. (B) Scatter plots showing coexpression patterns of Th17 and Treg lineage-specifying TFs (RORγt and Foxp3, respectively) under different input combinations: column 1, Th17-inducing conditions; column 2, IFNγ + TGFβ + IL-6 + IL-2; column 3, TGFβ + IL-4 + IL-2; column 4, iTreg-inducing conditions. (C) Mean levels of lineage-specific TFs measured after 7 and 14 d of culture under the indicated polarized (Left) and mixed (Right) conditions. Average normalized mean levels of two independent experiments are shown (SI Appendix, Extended Experimental Procedures). F, Foxp3; G, GATA3; R, RORγt; T, T-bet.
Organization of the CD4+ T-Cell Differentiation Space Reflects a Hierarchy of Cytokine Inputs.
Heterogeneity in T-cell responses was shown by a number of in vitro and in vivo studies to be averaged out for a large-enough initial cell population, such that the mean population response follows more deterministic rules (38, 43⇓–45). Therefore, we focus in the following analysis on the mean response of the cell population. To reduce dimensionality, we applied principal component analysis (PCA) to the measured mean protein levels. The first three principal components (PCs) account for 88% of the variability in the data and were used in our analysis. Organization of the measured 64 samples in the PC differentiation space is shown in Fig. 4 A and B. To allow for sample identification and comparison between data and model, we colored data points according to the similarity in the levels of the 10 proteins that we measured, as determined by sorting samples using the SPIN (sorting points into neighborhoods) algorithm (46) (SI Appendix, Fig. S4). Of note, the samples did not group into separate clusters in the PCA space but were spread over a wide region. Interestingly, the samples corresponding to the classical phenotypes (Th1, Th2, Th17, and iTreg) are found close to the edges of this region, whereas a continuum of phenotypes spans the range between them (Fig. 4 A–D).
Fig. 4.
Organization of cell states in the differentiation space reflects a hierarchy of input cytokines. PCA was performed on the mean expression levels of 10 measured proteins in all 64 experimental conditions. Each dot represents the response to 1 of 64 input conditions projected on the (A) PC1–PC2 plane and (B) PC1–PC3 plane. Samples are gradually colored (from dark blue to dark red; color bar on the right) according to the level of similarity in their protein expression pattern as determined using the SPIN algorithm (SI Appendix, Fig. S4). Samples are labeled according to input cytokine combinations using the one-symbol labels as specified in the table. Conditions that drive polarized Th lineages are marked: □, Th1; ♢, Th2; △, iTreg; ▿, Th17. (C) Hierarchy of input cytokines. The first PC divides the output differentiation space according to the presence of TGFβ as an input cytokine. The TGFβ+ subspace (PC1 < 0) is further divided by the third PC according the presence of IL-6. Samples that did not receive TGFβ as an input (PC1 > 0) were divided in a similar manner, this time according to the presence of IL-4. (D) A representation of the contribution of each measured output protein to cell differentiation. Each arrow represents the loading of a lineage-associated protein in the PC2–PC3 plane. The loading reflects the direction in which the change in the expression level of the corresponding protein is maximal. Loadings of most proteins point in the direction of the polarized state with which they are associated. (E) Relative expression levels of GATA3 (Left) and RORγt (Right) across the differentiation space indicated by color (color bar on the right). Conditions that drive polarized Th lineages are marked in C–E as described above.
The organization of cell phenotypes in the differentiation space (output) reflects the effect of the input cytokines in a hierarchical manner, with some inputs being more dominant than others. The first PC divides the space according to the presence of TGFβ as an input cytokine (Fig. 4C). Samples that received TGFβ as one of the input signals have negative PC1 values, whereas positive PC1 values were assigned to samples that were not cultured with TGFβ. The TGFβ+ subspace (PC1 < 0) is further divided by the third PC according the presence of IL-6. Samples that received both TGFβ and IL-6 as an input showed positive values on the PC3 axis, whereas samples cultured with TGFβ and without IL-6 had negative PC3 values. Samples that did not receive TGFβ as an input (PC1 > 0) were also divided by the third PC, this time according to the presence of IL-4 as an input (Fig. 4C). Similar hierarchy of input cytokines is also observed when samples are sorted according to similarity of the measured outputs using the SPIN algorithm (46) (SI Appendix, Fig. S4).
The organization in the differentiation space also reflects the level of correlation of each measured protein with the different cell states. We calculated for each measured output protein the direction in PC space in which its change is maximal (also termed the “loading” of each output). Most of the measured outputs have a positive loading along PC1 (SI Appendix, Fig. S5), reflecting a general suppressive effect of TGFβ, which lowers expression levels of most of the measured proteins, in accordance with previous reports (47) (SI Appendix, Fig. S5). When projected on the PC2–PC3 plane, the directions of output loadings are pointing toward the polarized phenotypes with which they are associated, thus reflecting the known association of lineage-specifying proteins with the polarized cell states (Fig. 4D). Notably, the four polarized phenotypes are located nearly orthogonal to each other in this projection, and the loadings of the associated proteins are also mostly orthogonal. Thus, GATA3, IL-4, IL-5, and IL-13 point in the general direction of the Th2-inducing condition IL-4 + IL-2 (diamond in Fig. 4D). IFNγ points in the opposite direction toward the Th1-inducing condition, IL-12. RORγt and IL-17 point in the direction of the Th17-inducing condition, TGFβ + IL-6, which is roughly perpendicular to the Th1–Th2 axis, and Foxp3 points in the opposite direction toward the Treg-inducing condition, TGFβ + IL-2 (Fig. 4D).
The mean expression levels of these lineage-specifying proteins change in a gradual way over the differentiation space, exhibiting high levels in the region that contains input conditions that support the corresponding lineage, which gradually drop away from this region (Fig. 4E).
A Hierarchical Additive Regression Model Depicts the Organization of CD4+ T-Cell Phenotypes in the Differentiation Space.
Based on the observed organization in the differentiation space, we aimed to develop a mathematical model that will relate the resulting phenotype, as assessed by the measured mean levels of lineage-specific proteins, to the input cytokine combinations. As a first approximation, we attempted to fit the data with a linear regression model:
Using our data, we solved the set of 64 linear equations per PC dimension and obtained seven parameters for each dimension (Ai,1,…, Ai,6, Bi) that best fit the data. Model accuracy was evaluated by the rms of the Euclidian distance between the model and the measured output (in the 3D PC space). Although the model fits the measured data significantly better than random shuffled data (P value < 0.001) (SI Appendix, Fig. S6), it does not capture the general organization of the measured phenotypes in the differentiation output space.
Based on the observed hierarchy of input cytokines (Fig. 4C), we hypothesized that the model can be improved by dividing the samples into groups based on shared inputs and fitting each group of conditions separately by a linear regression model. Therefore, we divided the 64 samples into four groups based on two input cytokines, Ia and Ib (a, b = 1,…, 6; a≠b): one group consists of all 16 samples that have both Ia and Ib as inputs, a second group consists of 16 samples that have neither of the two cytokines as an input, and two other groups have either Ia or Ib as an input. We then solved a linear model as described above for each group of 16 samples separately. In such a segmented linear regression model, there are different sets of parameters (Ai,1,…, Ai,6, Bi) that describe the linear input–output relations within each group. The parameters can differ between the groups, accounting for nonlinear transitions that are governed by the defined hierarchy (Ia, Ib).
We solved the hierarchical model for all of 15 possible combinations of two input cytokines (SI Appendix, Table S1). We also solved the model with single-input segmentation, when the data are divided into two groups based on the presence/absence of one of the input cytokines (SI Appendix, Table S1). Although the number of fitting parameters in the dual input segmentation is larger, not all two-cytokine combinations provided substantial improvement over the linear model or the single-input segmentation. For example, segmentation using the combination of [IL-2, IFNγ] (red in SI Appendix, Table S1) marginally improved model accuracy (Err = 1.129 compared with Err = 1.189 in the nonsegmented linear model and Err = 1.156 or 1.164 for single-input segmentation using IL-2 or IFNγ, respectively). Notably, three dual input segmentations, [TGFβ, IL-6], [TGFβ, IL-4], and [IL-4, IL-6], gave substantially better fit than the others (green in SI Appendix, Table S1). These three input cytokines are indeed those that had the strongest association with organization of cell phenotypes in the PC space (Fig. 4C). We note that our modeling approach does not allow us to favor one of these three best fit segmentations over the others. For consistency and because similar results were obtained using all three segmentations (SI Appendix, Fig. S7), we show in the following the results that were obtained using the [TGFβ, IL-6] segmentation (Err = 0.631; P value < 0.001; random shuffling test).
The hierarchical model captures the organization of the phenotypes in the differentiation space much better than the nonhierarchical one and reproduces the observed continuity of phenotypes ranging between the canonical Th states (Fig. 5 C and D and SI Appendix, Fig. S7). Model robustness was validated by a leave one out analysis (SI Appendix, Fig. S6) as well as by using only part of the samples to train the model and test its predictions on the rest of the samples that were not used for estimating model coefficients.
Fig. 5.
A segmented linear regression model describes CD4+ T-cell phenotypes obtained under diverse cytokine combinations. (A and B) Data: PCA was performed on the mean expression levels of 10 measured proteins under all 64 experimental conditions. Each dot represents the response to 1 of 64 input conditions projected on the (A) PC1–PC2 plane and (B) PC1–PC3 plane. Samples are colored as in Fig. 4. (C and D) Model: Results of a best fit hierarchical model using the [TGFβ, IL-6] segmentation. Results are projected on the (C) PC1–PC2 plane and (D) PC1–PC3 plane. Samples are colored as in A and B. (E) Graphical representations of the segmented linear model using vector addition in the PC differentiation space. Response to IL-4 (Left), IL-12 (Center), and their combination (Right). The response to both inputs is given by the sum of the vector responses to each input alone together with the “offset” response B, which describes contribution to the output that is not explicitly regulated by the input cytokines. (F) Comparing model predictions and measured protein levels. Parallel coordinates representation of the model solution (dots connected by dashed lines) and experimental measurements (diamonds connected by solid lines) of all 10 measured proteins in two representative input mixtures. Input cytokine combinations are marked using the one-symbol labels as specified in the table.
The model can be described graphically using vector algebra in the PC space. Fig. 5E shows the predicted response to three input conditions: IL-4, IL-12, and their combination IL-4 + IL-12. The response to the signal combination is given by the sum of responses to the individual signals (vector addition). Note that vector B is the same for all three conditions, because it describes the effective contribution of factors that do not depend explicitly on the inputs. This vector addition representation can also be extended to combinations of a larger number of input cytokines and serves as an intuitive graphical description of the input–output relations as revealed by our analysis.
In the above analysis, we solved the model in the 3D space defined by the three most dominant PCs. Because the model is solved for each dimension separately, we could also solve it in the 10D space of all measured output proteins. Solving in 10D space allows for a direct comparison between model predictions and the measured protein levels. Parallel coordinates representation (48) shows a very good agreement between model and data in the population-averaged levels of all 10 lineage-related proteins (Fig. 5F).
Together, these results show that the mean response of the cell population to combinations of cytokine inputs can be predicted based on the sum of the responses to individual inputs but only within the relevant segment. Different additive rules apply within different segments, reflecting the observed hierarchy of cytokine inputs.
The Segmented Linear Model Predicts Trajectories in the Differentiation Space in Response to New Input Mixtures.
The additive nature by which inputs are integrated according to the model suggests that new intermediate cell phenotypes can be achieved in a predictive way. To test this assumption, we attempted to populate the region in between Th1 and Th2 phenotypes in the differentiation space using varying concentrations of two input cytokines, IL-12 and IL-4. As we have previously shown (38), gradually changing concentrations of IL-12 and IL-4 result in mixed Th1–Th2 phenotypes. Fig. 6A shows the expression of IFNγ and IL-13 (Th1 and Th2 characteristic cytokines, respectively) for cells cultured under varying IL-12 and IL-4 concentrations. In the presence of high levels of IL-12 and no IL-4 (Th1 conditions), cells produce large amounts of IFNγ (Fig. 6A, column 1), whereas the production of IL-13 is induced in the presence of high levels of IL-4 and no IL-12 (Th2 conditions) (Fig. 6A, column 5). Under mixtures of two input cytokines, a mixed phenotype of cells coexpressing IFNγ and IL-13 is revealed (Fig. 6A, columns 2–4). The percentage of cells expressing these cytokines is changing in response to varying the concentrations of the input cytokines, gradually shifting between Th1 and Th2.
Fig. 6.
Model predicts differentiation outcome under new input conditions. (A) Flow cytometry measurements of IFNγ and IL-13 under mixtures of Th1- and Th2-inducing cytokines (IL-12 and IL-4, respectively) at varying concentrations. (B) Model predictions of cell state under eight mixtures of IL-4 and IL-12 provided at different concentrations. When maintaining a fixed concentration of one input cytokine while gradually changing the levels of the other, cell states move along a straight line in the PC space. (C) Experimental data: CD4+ T cells were cultured under eight mixtures of IL-4 and IL-12 at different concentrations, spanning the range between Th1 and Th2. Cell state, determined by mean expression levels of 10 measured proteins as above, is projected on the PC1–PC3 plane. Lines represent linear regression fits. In B and C, IL-12 input concentration is represented by the triangles on the left, ranging from dark to light red; IL-4 input concentration is represented by the triangles on the right, ranging from dark to light blue (as indicated by the color bars on the right); all 64 conditions of Fig. 4 are shown as gray dots as a reference. (D) The 3D representation of the experimental data of C from three different points of view of the PC space. All eight points fall close to one plane in the PC space as predicted by the segmented linear model. (E) Model predictions of protein expression patterns for eight combinations of IL-12 and IL-4 at varied concentrations as in C. Input cytokine combinations are indicated by the color bars below using the same color scheme as in B and C. (F) Measured mean protein expression levels of cells cultured under the same conditions as in C.
As we measure the expression levels of a larger number of proteins, we can gain a broader view of how cell state shifts in response to changes in input cytokine levels. Thus, we tested the projection of the mixed phenotypes in the higher-dimensional PC differentiation space, which describes the population-averaged levels of all 10 measured proteins. According to the linearity of the model, changing the concentration of one input while maintaining the concentrations of other inputs fixed should result in phenotypes that are positioned along a straight line in the PC differentiation space. Thus, the model predicted that increasing IL-4 levels from zero to high concentrations while maintaining IL-12 at a fixed high concentration will generate intermediate phenotypes along a line that starts at the Th1 region and approaches a region in between Th1 and Th2 in the PC space (Fig. 6B). Accordingly, when reducing IL-12 levels from high to zero while maintaining high concentration of IL-4, the model predicts intermediate phenotypes along a line at a different angle, ranging from the “mixed” Th1–Th2 region to the Th2 region (Fig. 6B).
To test these predictions experimentally, we cultured cells under mixed Th1–Th2 input condition. Four culture conditions had a fixed high level of IL-12 and increasing levels of IL-4. Additionally, four samples had a fixed high level of IL-4 and decreasing levels of IL-12 (down to zero). The same 10 output proteins were measured as in the previous experiments. When projected onto the PC1–PC3 plane, the resulting measured phenotypes closely follow two straight lines (Fig. 6C) as predicted by the segmented linear model (Fig. 6B). Furthermore, the angles of these two lines in the PC space are in agreement with those predicted by the model. Of note, the model predicts that all of these eight points will fall on a single plane in PC space, which is defined by two linear lines. Indeed, we find that the measured states all fall very close to one plane in the PC space (Fig. 6D). The fact that all eight points fall on one plane presents a tight restriction of the possible degrees of freedom, which is predicted by the model and confirmed by the experimental data.
We also used the model to calculate the predicted mean levels of all 10 measured proteins under the eight Th1–Th2 mixtures. Model predictions of protein levels (Fig. 6E) are in general good agreement with their experimentally measured values (Fig. 6F). This analysis provides additional insights regarding the input–output relations of the system as described by our model. Levels of the Th2-associated proteins, such as GATA3 and IL-13, gradually increase along the first line (Th1 to mix; fixed IL-12, increasing IL-4) and remain almost constant along the second line (mix to Th2; decreasing IL-12, fixed IL-4). Levels of IFNγ, which is Th1-associated, gradually decrease along the first line (Th1 to mix) and remain almost constant along the second line (mix to Th2). These trends are predicted by the model and were observed in the experimental data. Notably, the expression of IL-10 shows a nonmonotonic behavior—its level under mixed conditions is higher than in either Th1 or Th2, which was predicted by the model and verified experimentally. This trend was also observed in the initial experiments of 64 input conditions. Thus, the model provides quantitative predictions of how the expression levels of lineage-associated proteins change in response to changes in inputs concentration, even for proteins that are not conventionally associated with the specific input cytokines.
Discussion
In this study, we asked how naïve CD4+ T cells combine a number of cytokine inputs to determine their differentiation state. We studied this question by systematically exposing naïve CD4+ T cells to varied cytokine combinations and measuring their resulting response in terms of changes in expression of key genes. We found that, under most conditions, rather than committing to a polarized fate, CD4+ T cells adopted varied mixed and intermediate phenotypes and expressed diverse combinations of lineage-specific cytokines and TFs. Moreover, we reveal that the expression level of lineage-specific proteins in a population of cells is heterogeneous and spans a range of intermediate levels that are not limited to extreme values (either high or low). Cell populations with mixed phenotypes, displaying the characteristics of different lineages simultaneously, may increase the flexibility of the immune response in complex cytokine environments. Although previous studies have reported coexpression of lineage-specifying TFs and lineage-characteristic cytokines (36, 39, 49⇓⇓⇓–53), the rules by which combinations of input signals lead to various patterns of coexpression remained elusive.
Combined with mathematical modeling, our data allowed us to determine the rules that the cells follow when responding simultaneously to a number of input signals. We found that cell populations respond to cytokine combinations in a way that is additive in nature, such that the mean response to a combination of cytokines equals the sum of the responses to the individual cytokines. Additive signal integration was observed in other studies of cellular responses (45, 54⇓⇓⇓–58). We show that, in the case of CD4+ T-cell differentiation, however, there is a hierarchy that modulates the summation of signals: some inputs are more dominant, and if present, they modify the coefficients of the additive model. For example, in the absence of TGFβ, there is one set of coefficients for the summation of the other inputs; in the presence of TGFβ, there is a different set of coefficients, but the response still remains additive in nature.
Our modeling shows that dividing the conditions into four groups according to the presence or absence of two dominant inputs, while maintaining additive responses within each group of inputs, captures the behavior of the measured cell phenotypes. The segmented linear regression model provides a general framework for evaluating the level of nonadditivity and identifying the inputs that influence it the most. Additional studies may explore the biological mechanisms that lead to the two-level hierarchy and reveal mechanistic reasons for favoring a specific segmentation over the others.
The additive nature of the response is surprising, because gene regulation during Th differentiation involves nonlinear interactions with many positive and negative feedback loops. The relatively simple rules that we propose here suggest that the underlying gene regulatory networks have a limited level of functional complexity in their input–output relations, which is much lower than their observed topological complexity. Bridging the gap between the growing understanding of complex and nonlinear regulatory networks that coordinate cell differentiation and the observed simpler additive rules for cellular responses remains a main challenge for systems biology and systems immunology.
Previous models described differentiation of Th cells toward a specific fate (55⇓–57) or focused on key molecular components (38, 58⇓–60). Other models describe Th differentiation at a higher level of detail by modeling regulatory networks of increasing complexity (61⇓⇓–64). These different modeling approaches promoted our quantitative understanding of Th cell differentiation. However, because the majority of experimental data are qualitative, the scarcity of quantitative data for most interactions remains a main challenge for modeling Th differentiation and other complex biological systems based on underlying molecular interactions. Our black box approach provides an alternative modeling scheme, because it describes input–output relations based on extensive experimental data without specifying the response of every node in the regulatory network. This approach allowed us to study the logic of operation of the differentiation process and reveal rules that describe the mean cellular response. Furthermore, our model could successfully predict cell phenotypes measured in additional experiments, which were not used for its construction. This approach can be further combined with studies that map signaling responses to physiological outcomes (65, 66), leading to improved models for predicting cellular responses to input combinations.
There are several samples that reside near the origin of the PCA differentiation space (Fig. 4 A and B). This region is indicating intermediate levels of expression of most of the measured proteins. One explanation may be that the input cytokines in these samples are antagonistic and cancel the effect of each other. It is also possible that some samples that reside in this region express other proteins that we did not measure in this study. For example, it was shown that a combination of TGFβ + IL-4 drives differentiation of Th9 cells (67, 68), which can be revealed by also measuring IL-9 as an output. New technologies, such as mass cytometry (69) and single-cell RNA-seq (70, 71), will allow us to extend these studies to higher dimensions. The hierarchical additive model can be used to describe systems with higher numbers of outputs, because the model is solved independently for each dimension.
We find that coexpression of lineage-specifying TFs is stable for at least 14 d of culture under mixed input conditions. Coexpression patterns did not become more polarized during a prolonged culture and may even show a more mixed phenotype in some cases (Fig. 3C). These results suggest that the mixed states do not represent differentiation intermediates that eventually converge into one of the polarized states. This notion is supported by previous data on mixed Th1–Th2 cells. We and others observed stability of a mixed Th1–Th2 phenotype for the duration of weeks (in vitro) (38) and months (in vivo) (39). Devising adequate in vivo models to test for stability and immune function of other mixed Th phenotypes is challenging but can be insightful for better understanding of their physiological significance.
The ability of our model to predict cellular responses to complex cytokine combinations can provide a framework for understanding cellular behavior and predicting responses to new cytokine combinations and doses. The use of such quantitative approaches can assist in the development of immunotherapy strategies that use ex vivo activation and differentiation of CD4+ T cells. Potential applications include ex vivo generation of T cells with a desired phenotype for adoptive cell therapy and the use of cytokine combinations for rational design of new therapeutic approaches.
Experimental Procedures
Mice.
C57BL/6 female mice where purchased from Harlan Laboratories. Mice were housed under specific pathogen-free conditions at the animal facility of the Weizmann Institute and used at 6–8 wk of age. All animal experiments were performed under protocols approved by the Animal Care and Use Committee of the Weizmann Institute.
Cell Culture.
Naïve CD4+ T cells were isolated from spleens of mice by magnetic microbeads (Miltenyi Biotec) and cultured with plate-bound anti-CD3 and anti-CD28 at 1 × 106 cells per 1 mL concentration. Similar results were obtained by FACS sorting CD4+CD25−CD62LhighCD44low naïve T cells (SI Appendix, Fig. S8). Cells were supplemented with different combinations of cytokines used in the following concentrations unless stated otherwise: IL-12 (10 ng/mL), IFNγ (10 ng/mL), IL-4 (20 ng/mL), TGFβ (10 ng/mL), IL-6 (20 ng/mL), and IL-2 (5 ng/mL; all from R&D). Full materials are given in the extended experimental procedures in SI Appendix.
Cell Proliferation Assay.
Naïve CD4+ T cells were isolated from spleens of mice by magnetic microbeads (Miltenyi Biotec) and labeled with 2 µM carboxyfluorescein diacetate succinimidyl ester (CFSE; Invitrogen) for 10 min at 37 °C. After incubation, cells were washed with PBS, resuspended in RPMI, and cultured with plate-bound anti-CD3 and -CD28 at a concentration of 1 × 106 cells per 1 mL. Cells were supplemented with the appropriate cytokines and antibodies as indicated. Cell proliferation was assessed by CFSE dilution after 72 h of culture using flow cytometry and further analyzed with FlowJo software (TreeStar).
Antibodies, Surface, and Intracellular Staining.
For intracellular staining, cells were restimulated 7 d after their initial stimulation with plate-bound anti-CD3 and stained with anti–IL-4, –IL-5, –IL-13, –IL-17, -IFNγ, –IL-10, -GATA3, –T-bet, -Foxp3, -RORγt, and -CD4. Samples were measured using flow cytometry (BD). Full materials are given in the extended experimental procedures in SI Appendix.
Mass Cytometry (CyTOF).
After a 7-d culture, cells were restimulated and stained for various surface and intercellular proteins as specified in SI Appendix, Table S2. Cell barcoding and preparation for mass cytometry were conducted as previously described (72, 73). Samples were measured using mass cytometry (CyTOF Helios mass cytometer; Fluidigm). See extended experimental procedures in SI Appendix for full details.
RNA Preparation, Library Construction, and RNA Sequencing.
After culture with specified cytokine mixtures, 1 × 105 cells were lysed in lysis binding buffer, and polyadenylated RNA was purified using oligo dT magnetic beads (Ambion). Libraries were constructed according to the MARS-seq (massively parallel RNA single-cell sequencing) protocol (71) adjusted to produce low-input RNA-seq libraries for the Th cell populations studied in this work. Primers used in MARS-seq are summarized in SI Appendix, Table S3. Full materials are given in extended experimental procedures in SI Appendix. Each sample was sequence to the depth of 5 million reads using an Illumina HiSeq 1500.
Data Normalization and Analysis.
Mean fluorescence levels (MFI) were standardized by subtracting for each protein its mean level averaged across all tested input conditions and dividing by the SD of the MFI values across all input conditions. Arithmetic mean was used. Using a geometric mean gives very similar results (correlation of R = 0.98). The median values of four independent experiments were used in the analysis (SI Appendix, Fig. S8 shows comparison between independent experiments). Expression levels of mRNA were standardized in a similar manner. In the analysis presented in Fig. 2 G–J, expression levels were normalized by the level of the lowest measured sample (in a log scale). In the 14-d experiments (Fig. 3C), data were normalized, such that the level of each TF in the condition that supports its associated lineage equals one, and its average level under three other polarizing conditions is zero. Data analysis was conducted using Matlab (MathWorks), Easyflow (a Matlab package that was developed in our laboratory to analyze flow cytometry data), and the SPIN algorithm (46).
Mathematical Modeling and Statistical Analysis.
Linear model.
We solved a linear regression model as described by Eq. 1. Because the set of equations is underdetermined (64 equations with seven parameters for each PC dimension), we found the least squares solution to the system of equations using Matlab. We implemented the solution using the \ (“mldivide”) function in Matlab, which uses QR (orthogonal–triangular) decomposition for this purpose. Model accuracy was evaluated by the rms of the Euclidian distance between the measured output vector and the model estimations. Model results were evaluated using a leave one out analysis as well as the data shuffling test (SI Appendix, Fig. S6).
Segmented linear model.
We divided all samples into four groups of 16 samples each, based on a hierarchy of two input cytokines, Ia and Ib (a,b = 1,…, 6; a≠b) as described. We then solved a linear model as described above for each group of 16 samples separately. Within each hierarchy, we have five unknown coefficients (corresponding to four cytokines that do not define the hierarchies and the “bias” coefficient Bi) and 16 equations. Thus, four systems of equations are underdetermined and were solved as described above.
Varied cytokine concentrations (Fig. 6) were introduced into the model by varying values in the I input matrix using the following mapping: 540 ng/mL = 2.5, 90 ng/mL = 1.5, 15 ng/mL = 1, and 2.5 ng/mL = 0.5.
Acknowledgments
We thank members of the laboratory of N.F. for helpful comments and discussions. This research was supported by Israel Science Foundation Grant 1184/15 and the David and Fela Shapell Family Foundation Israel National Center for Personalized Medicine Fund for Preclinical Studies.
Footnotes
↵1I.E.-M. and J.R. contributed equally to this work.
- ↵2To whom correspondence should be addressed. Email: Nir.Friedman{at}weizmann.ac.il.
Author contributions: I.E.-M., J.R., and N.F. designed research; I.E.-M., J.R., I.Z., D.L.-A., and S.R.-Z. performed research; I.E.-M., J.R., and N.F. analyzed data; and I.E.-M., J.R., and N.F. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission. F.S. is a guest editor invited by the Editorial Board.
Data deposition: The flow cytometry data in this study are available upon request.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1615590114/-/DCSupplemental.
References
Hierarchical summation of signals by CD4 T cells
Sign up for Article Alerts
Article Classifications
- Biological Sciences
- Systems Biology
Jump to section
You May Also be Interested in
Astronomers thought they’d finally figured out where gold and other heavy elements in the universe came from. In light of recent results, they’re not so sure.
Image credit: Science Source/Tom McHugh.
Researchers are still trying to understand what causes this strong correlation between neural and social networks.
Image credit: Shutterstock/Yeongsik Im.
Amanda Rodewald, Ivan Rudik, and Catherine Kling talk about the hazards of ozone pollution to birds.
CRISPR-Cas9 gene editing can improve the effectiveness of spermatogonial stem cell transplantation in mice and livestock, a study finds.
Image credit: Jon M. Oatley.