Research ArticleImmunology

Cell-to-Cell Variability Analysis Dissects the Plasticity of Signaling of Common γ Chain Cytokines in T Cells

See allHide authors and affiliations

Science Signaling  12 Mar 2013:
Vol. 6, Issue 266, pp. ra17
DOI: 10.1126/scisignal.2003240


Natural variability in the abundance of signaling regulators can lead to divergence in cell fate, even within genetically identical cells that share a common differentiation state. We introduce cell-to-cell variability analysis (CCVA), an experimental and computational methodology that quantifies the correlation between variability in signaling regulator abundance and variation in the sensitivity of cells to stimuli. With CCVA, we investigated the unexpected effects of the interleukin 2 (IL-2) receptor α chain (IL-2Rα) on the sensitivity of primary mouse T lymphocytes to cytokines that signal through receptors that have the common γ chain (γc). Our work showed that increased IL-2Rα abundance decreased the concentration of IL-2 required for a half-maximal activation (EC50) of the downstream effector signal transducer and activator of transcription 5 (STAT5), but reduced the responsiveness to IL-7 or IL-15, without affecting the EC50 values of other γc cytokines. To investigate the mechanism of the effect of IL-2Rα on γc cytokine signaling, we introduced a Bayesian-inference computational framework that models the formation of receptor signaling complexes with data from previous biophysical measurements. With this framework, we found that a model in which IL-2Rα drives γc depletion through the assembly of functional IL-2R complexes was consistent with both the CCVA data and experimental measurements. The combination of CCVA and computational modeling produced quantitative understanding of the crosstalk between γc cytokine receptor signaling in T lymphocytes.


Quantifying the effect of protein abundance on cellular function has attracted considerable attention in recent years (14). To perform this analysis in bacteria, researchers change incrementally the abundance of a chosen protein and then measure the functional consequences (59). However, this approach is more cumbersome in primary mammalian cells, such that protein function has been principally studied in an all-or-nothing fashion with genetic mutants or RNA interference. As an alternative, we propose that the inherent natural variability in protein abundance, as observed within populations of genetically identical mammalian cells (1015), can be used to dissect the quantitative regulation of signal transduction.

To assess the phenotypic variability of populations of isogenic cells, researchers can quantify the variability in the abundance of mRNA or protein with individual cell resolution with single-cell quantitative reverse transcription–polymerase chain reaction assays (16) or flow cytometry (12, 13). Notably, studies using these techniques have demonstrated large heterogeneity in the abundance of signaling components (including receptors, kinases, adaptors, phosphatases, and cytokines), with typical coefficients of variation (CVs) for the lognormal distribution of mRNA and protein amounts larger than 0.5 within activated T cell clones, for example (16, 17). Thus, in such distributions, 15% of cells will have protein abundances that deviate from the median by more than twofold. Even larger variability was uncovered in the case of the interleukin-2 (IL-2) receptor α chain (IL-2Rα), with CVs of up to 3.0 in populations of genetically identical transgenic T cells activated in vitro (17). In these cells, 15% of the population of cells has an IL-2Rα abundance that deviates from the median by more than 10-fold in either direction (17). In settings of infection, this variability in the abundance of IL-2Rα on T cells correlates with a split between short-lived effector and memory precursor fates (18). Examples in which a continuous spectrum of surface protein abundance maps onto discrete differentiation paths have been reported in other biological systems as well (12, 19). These observations raise the question of how variability in protein abundance affects signaling thresholds and ultimately cell differentiation decisions.

Here, we present an experimental methodology to quantitatively correlate such variability in protein abundance with variable regulatory function. As part of this methodology, we present a software package, ScatterSlice, which quantitatively correlates the functional heterogeneity of cellular signaling thresholds with variation in protein abundance. We applied this methodology to investigate the signaling responses of T cells to the common γ chain (γc) cytokines—IL-2, IL-4, IL-7, IL-9, IL-15, and IL-21—whose receptors share the γc chain (also known as CD132) (20). The receptors for IL-2 and IL-15 also share their β chain (CD122). Each cytokine receptor has an α chain that imparts specificity through enhanced affinity for its respective cytokine. For example, the IL-2R is formed through the assembly of the γc and β chains with a specific α chain, IL-2Rα (also known as CD25), whereas the IL-7R is composed of the γc and its specific α chain, IL-7Rα (also known as CD127). Although it has been established that the α chain of each cytokine receptor is required for signaling (21), a quantitative understanding of the effects of sharing a common receptor chain on the ability of T cells to trigger downstream phosphorylation in response to γc cytokines remains elusive.

Specifically, we investigated whether variation in the abundance of IL-2Rα influenced the sensitivity of T cells to IL-2 and other γc family cytokines. In agreement with previous results (17), we report that a large variation in IL-2Rα abundance correlated with much of the variability in the sensitivity of cells to IL-2. More surprisingly, we found that a higher abundance of IL-2Rα correlated with inhibited sensitivity of the cells to IL-7. We ruled out decreases in IL-7Rα abundance at the cell surface or direct inhibition of IL-7Rα signaling as a mechanism to explain the dependency of IL-7 sensitivity on IL-2Rα. Instead, a simple thermodynamic model suggested that the large abundance of IL-2Rα could inhibit IL-7 signaling indirectly by sequestering γc chains into preformed heterotrimeric IL-2 receptors in the absence of IL-2. Extending our experimental system to other cytokines of the γc family, we determined that the sensitivity of cells to IL-15 was slightly dampened by increased IL-2Rα abundance, whereas sensitivities to IL-4 and IL-21 were unaffected. We accordingly expanded our computational model to the entire γc family to investigate how varied abundance of IL-2Rα inhibited cellular responsiveness to IL-7 and IL-15, but not IL-4 and IL-21. This analysis indicated that the affinity between the α chain of each cytokine receptor (for example, IL-21Rα) and the γc chain was the critical factor determining whether increased IL-2Rα abundance could inhibit sensitivity to that cytokine. These studies illustrate the usefulness of combining high-resolution data from cell-to-cell variability analysis (CCVA) with computational modeling to gain a quantitative understanding of cellular signaling.


Phenotypic variability between T cell populations leads to differences in IL-2 EC50 values

We sampled isogenic populations of primary 5C.C7 T cell receptor (TCR)–transgenic RAG−/− IL-2−/− mouse T cells at different times after activation with congenic B10.A CD3ε−/− splenocytes pulsed with 1 μM K5 antigenic peptide, and found that a given concentration of IL-2 stimulated greatly varied signaling responses, as reflected in the different degrees of phosphorylation of signal transducer and activator of transcription 5 (STAT5) (Fig. 1A). We then analyzed the abundance of phosphorylated STAT5 (pSTAT5) in T cells in response to a complete range of concentrations of IL-2 after a 10-min exposure. In one culture, we estimated the typical concentration (EC50) of IL-2 that yielded a half-maximal pSTAT5 response to be 10 pM (Fig. 1B, green line), which is consistent with the estimated dissociation constant (KD) for the “high-affinity IL-2 receptor” (22) as well as with functional readouts, such as the proliferation of CTLL-2 cells (23). However, in other cultures, we observed shifts in the EC50 for IL-2 by 10-fold—up to 100 pM or down to 1 pM (Fig. 1B). Flow cytometry measurements of IL-2Rα abundance revealed a close inverse correlation between the IL-2 EC50 and the median abundance of IL-2Rα (Fig. 1C). Cells from a culture with a higher EC50 had a lower median abundance of this receptor subunit, whereas cells from a culture with a lower EC50 had a greater median abundance. Overall, the median abundance of IL-2Rα per cell correlated well with the IL-2 EC50 (R2 = 0.87; Fig. 1C).

Fig. 1

IL-2 EC50 values are inversely correlated with IL-2Rα abundance. In vitro cultures of primary mouse CD4+ T cells were examined at day 3 (pink), day 5 (green), day 6 (orange), or day 30 (blue) after initial in vitro activation with peptide–major histocompatibility complex antigen and proliferation in IL-2–supplemented medium. Cultures were selected from 10 experiments to show the range of responses. (A) Ten minutes after stimulation with 10 pM IL-2 (colored histograms), T cells display variability in their distribution of pSTAT5. These distributions were normalized between 2 and 98% of the unstimulated and maximal responses of each culture, respectively. The gray histogram depicts the abundance of pSTAT5 in unstimulated (control) cells. (B) Hill equations (curves) fit to normalized dose responses (filled circles) display large variability in EC50 values (crosses), the characteristic concentration of IL-2 for which 50% of the maximal pSTAT5 response is reached. (C) The EC50 for the response of cells to IL-2 varies over two orders of magnitude and is inversely correlated with the abundance of IL-2Rα. Black squares represent additional cultures of T cells not depicted in other panels. Error bars show the error in EC50 fitting. (D) Histograms of calibrated IL-2Rα abundance for each of the cell cultures under consideration. (E) Scheme for the dual labeling of IL-2Rα with different antibodies (F) Flow cytometric analysis of cells singly labeled with PE-Cy7–conjugated (red) or fluorescein isothiocyanate (FITC)–conjugated (blue) antibodies for IL-2Rα demonstrates a CV of 1.5, reflecting a wide range of per-cell IL-2Rα abundances within the population. In cells incubated with both antibodies simultaneously (purple), however, the ratio of the two stains has a CV of only 0.2, demonstrating that independent flow cytometric measurements of IL-2Rα abundance concur closely.

Cell-to-cell variation in IL-2Rα abundance leads to intrapopulation variation in EC50

Cells within a given population had differences in IL-2Rα abundance of the same magnitude as those in median IL-2Rα abundance between different cell populations (Fig. 1D). Because the variability in IL-2 EC50 between T cell cultures correlated well with the variability in median IL-2Rα abundance, we conjectured that the same correlation should exist when examining individual cells within a population.

This cell-by-cell assessment required that IL-2Rα abundance be measured with precision by flow cytometry. Dual labeling of IL-2Rα with two different antibodies (Fig. 1E) revealed a broad range of abundance observed within each fluorescence channel but minimal variation in the ratio of the two fluorescence measurements (Fig. 1F). Thus, even with the combined noise from staining with two antibodies and measurement across two channels, IL-2Rα abundance, which varied by up to 105-fold, was determined with a less than 1.5-fold error. Cell sorting experiments further confirmed that distinct subpopulations were accurately defined by their IL-2Rα abundance (fig. S1).

To quantify systematically the correlation between cytokine EC50 and receptor abundance, we developed a three-stage methodology, named CCVA (Fig. 2A). First, cells were treated with a range of concentrations of stimulus. Second, protein abundance and downstream responses were measured concomitantly with single-cell resolution by flow cytometry. Finally, the flow cytometry data were analyzed with a software package named ScatterSlice (Fig. 2, B to E). Conventional software programs for the analysis of flow cytometry data deliver snapshots of the cellular response for a given concentration of stimulus, with the possibility of manually gating for subpopulations based on protein abundance. In contrast, our software was specially designed to automatically parse the heterogeneous population into subpopulations as defined by protein abundance (Fig. 2B), quantify the extent of phosphorylation of downstream signaling molecules in each subpopulation (Fig. 2C), and then determine the sensitivity of each subpopulation at all concentrations of a given stimulus to derive its EC50 within each subpopulation as a measurement of its stimulus sensitivity (Fig. 2D). Overall, CCVA delivers a complete map of the relationship between protein abundance and response sensitivity as quantified by the EC50 for the stimulus (Fig. 2E).

Fig. 2

CCVA methodology. (A) Outline of the overall methodology. Cells are first treated with a range of concentrations (dose response) of the stimulus, which can be a cytokine, drug, or any treatment that results in changes in cells that are measurable by flow cytometry. After a short stimulation, cells are processed for analysis by flow cytometry. After acquisition, the flow cytometry data are processed through ScatterSlice, a graphical user interface for CCVA. (B to E) Overview of the internal processing of data with ScatterSlice. (B) Binning the population for different amounts of IL-2Rα. The colors of the dots are carried over to (D) and (E). (C) The median value of the response, in this case the abundance of pSTAT5, is calculated within each IL-2Rα bin for each concentration of stimulus. Open histograms, no cytokine control; red histograms, 10 pM IL-2. (D) Within each bin, across all concentrations of stimulus (IL-2), a curve with the parameters of EC50 and amplitude (solid lines) is fit across all IL-2 concentrations for which a sufficient number of cells exist within a bin to accurately estimate the amount of pSTAT5 (circles and dotted lines). (E) The output of the curve fitting is the relationship between IL-2Rα abundance and EC50 (colored circles). The original histogram (gray) is shown for reference. Error bars show uncertainty in EC50 fitting.

We applied CCVA to quantify the effect of variation in IL-2Rα abundance on the IL-2 EC50 within one population of T lymphocytes. We found that as receptor abundance increased from 102 to 107 proteins per cell, the EC50 of the IL-2 response decreased from 100 to 1 pM (Fig. 2E). Extending the same analysis to different T cell populations (Fig. 1) revealed that despite large differences in median IL-2Rα abundance and corresponding shifts in the overall IL-2 EC50 from population to population (Fig. 1C), an identical inverse correlation between IL-2Rα abundance and EC50 was observed (fig. S2).

CCVA uncovers an increased IL-7 EC50 that correlates with increased IL-2Rα abundance

In effector T cells, IL-2Rα (also known as CD25) is commonly used as a marker of activation, identifying an active effector phenotype between the naïve and memory stages of differentiation (24). Furthermore, effector T cells with a lower abundance of IL-2Rα preferentially commit to a memory phenotype, whereas those with a higher abundance of IL-2Rα do not (18). Because both naïve and memory phenotypes are classically associated with enhanced sensitivity to IL-7 (25), we used CCVA to quantify the relationship between IL-2Rα abundance and the IL-7 EC50. This analysis demonstrated a positive correlation between IL-2Rα abundance and IL-7 EC50, with increased IL-2Rα abundance (from 102 to 107 receptors per cell) scaling with increased EC50 for IL-7 (from 3 to 60 pM) (Fig. 3A).

Fig. 3

Positive correlation between the EC50 for the IL-7 response and IL-2Rα abundance: CCVA and modeling. (A) IL-7 EC50 (red) and IL-2 EC50 (blue) as functions of IL-2Rα abundance, as measured by CCVA. Points show geometric means ± SD of EC50 values from CCVA of two independent dose responses with two replicates. These data are representative of 10 independent experiments. (B) Assembled cytokine-receptor complexes for IL-2 and IL-7. (C) Schematic of the selected model. Closed circles, single components; open circles, multimolecular complexes; colored circles, complete IL-2 (blue) and IL-7 (red) signaling complexes, consisting of IL-2–IL-2R and IL-7–IL-7R, respectively. A given complex sits on the line joining the two reactants involved in the corresponding binding reaction. (D) Comparison between the data and the fit for the model represented in (C) for different constraints on the parameters. Left panel, model output with measured parameters Kexp; middle panel, fit with constraints; right panel, unconstrained fit. (E) Modeled number of preformed IL-2Rs and IL-7Rs (in the absence of extracellular cytokine) as a function of IL-2Rα abundance for the model represented in (C) with the parameters obtained from the constrained fit (D, middle panel). (F) Cells were sorted into two groups, the 35% of the population that had the greatest abundance of IL-2Rα and the 35% of the population that had the least abundance of IL-2Rα, after which they were subjected to immunoprecipitation (IP) reactions with antibody against the γc chain and Western blotting analysis for IL-7Rα. These data show diminished IL-7Rα abundance in complex with γc in cells within the population that had the greatest abundance of IL-2Rα. The Western blot shown is representative of two independent experiments.

A decrease in IL-7R abundance or signaling does not account for the observed increase in IL-7 EC50

The effect of IL-2Rα abundance on IL-7 EC50 could be hypothesized to stem from the reported inverse correlation between IL-2Rα and IL-7Rα abundances during the course of differentiation from naïve to effector to memory phenotypes (25, 26). However, in our system, flow cytometric analysis demonstrated that the abundances of IL-2Rα and IL-7Rα at the cell surface were not inversely correlated (fig. S3A). We confirmed this result by sorting cells according to IL-2Rα abundance and demonstrating by Western blotting a lack of correlation with IL-7Rα abundance (fig. S3B). More critically, our investigations of EC50 at the population level revealed that although the IL-7 EC50 varied from 3 to 60 pM between different populations of cells, this variation did not correlate with IL-7Rα abundance with any statistical significance (fig. S3C). This lack of correlation reflects a lack of uniformity in EC50 despite fairly uniform abundance of IL-7Rα across effector cell cultures, suggesting that factors other than IL-7Rα abundance regulated the IL-7 EC50 in effector cells.

As an alternative mechanism, we tested whether IL-2Rα abundance correlated with a direct inhibitory modification of the IL-7Rα chain, such that cells with equal abundance of IL-7Rα could nonetheless differ in the number of functional IL-7Rs. If this were the case, all cytokines signaling through IL-7Rα would have a reduced EC50 that correlated with the higher abundance of IL-2Rα. To test this possibility, we applied CCVA to characterize whether the EC50 of T cells for thymic stromal lymphopoietin (TSLP), a cytokine whose receptor also contains IL-7Rα (fig. S4A), also correlated with IL-2Rα abundance. Our results showed no correlation between TSLP EC50 and IL-2Rα abundance (fig. S4B). Therefore, the increase in IL-7 EC50 could not be attributed to direct inhibition of IL-7Rα signaling.

An equilibrium binding model of γc chain competition quantitatively accounts for the effects of IL-2Rα on IL-7 EC50

We then focused on the fact that IL-2R and IL-7R share the γc chain (Fig. 3B) (21, 27). Specifically, we hypothesized that the increased abundance of IL-2Rα might drive the preformation of complete heterotrimeric IL-2 receptors composed of IL-2Rα, IL-2Rβ, and γc chains. This sequestration of the γc chain in preformed IL-2R complexes might lead to the increase in IL-7Rα EC50 observed in cells with increased IL-2Rα abundance (Fig. 3A).

Experimental evidence for the existence of complete heterotrimeric IL-2R complexes on the surfaces of live cells in the absence of cytokine remains equivocal. Indeed, the crystal structure of IL-2 in complex with the soluble extracellular portions of the full receptor indicates no contact between IL-2Rα and either the β or the γc chain (28). Moreover, isolated extracellular domains of IL-2Rα and γc have no detectable binding affinity in solution (29). On the other hand, the extracellular domains of IL-2Rα and the β chain have a binding affinity of 280 nM, as measured in solution (29). Furthermore, all pairwise interactions between IL-2Rα, the β chain, and the γc chain have been demonstrated on the surface of live cells in the absence of IL-2 (20, 30). It is therefore possible that in the absence of IL-2, an increased abundance of IL-2Rα could drive the preformation of full heterotrimeric IL-2R complexes. We therefore investigated whether this hypothesis could realistically and quantitatively account for the correlation between IL-2Rα abundance and IL-7 EC50 observed in our CCVA experiments by constructing an equilibrium binding model of IL-2R and IL-7R assembly and signaling. We first computed the EC50 of the cellular response to IL-2 and IL-7 from affinity constants (see the Supplementary Materials for more details) and quantitatively compared these results to values measured by CCVA. We then identified the simplest equilibrium binding model with realistic affinity constants that could account for the correlation between IL-2 and IL-7 EC50 values with the abundance of IL-2Rα as uncovered by CCVA.

Estimating the cytokine EC50 with an equilibrium binding model for the receptor chains is appropriate because (i) the system reaches a steady state within 10 min of exposure to cytokines (fig. S5); (ii) the abundance of pSTAT5 correlates linearly with receptor occupancy (fig. S6); and (iii) depletion of cytokines is negligible (fig. S7), such that cytokine concentrations can be considered to be constant during our experimental time scale. In this framework, a model for cytokine EC50 is determined by the set of possible receptor complexes formed or, equivalently, by the set of affinity constants associated with their formation (see the Supplementary Materials for more details). For this system, the abundances of all of the cytokine receptor chains can be measured (fig. S8, and see Materials and Methods), and all of the potentially required affinity constants have been experimentally determined, directly or indirectly (fig. S9) (20, 29, 3136). Thus, one should be able to create a model reproducing the data directly from these experimentally determined affinity constants. On the basis of the potential interactions between subunits of IL-2R and IL-7R, there exist 553 relevant models that can be enumerated systematically (figs. S10 and S11, and see the Supplementary Materials). However, none of the possible models quantitatively reproduces the positive correlation between the IL-7 EC50 and IL-2Rα abundance when affinity constants are fixed to their experimentally determined values (Fig. 3D, “Model with measured affinity constants”). (The best among these poorly performing models is shown in fig. S12.) Alternatively, when the affinity constants are unconstrained and free to vary, very good agreement between data and model output can be obtained for a large number of models (Fig. 3D, “Model with free affinity constants”) but at the expense of fitted affinity constants that greatly deviate from their measured values (fig. S13).

This discrepancy can be resolved by acknowledging the presence of some intrinsic uncertainties in the experimental determination of affinity constants (see fig. S9) and in the translation of in vitro measurements to live cells. To do so, we introduced a Bayesian strategy to model the CCVA-measured EC50 values while relying on previous measurements of affinity constants to constrain the fitting of model parameters (see the Supplementary Materials for more details). This procedure enabled us to identify a set of models that achieved good agreement with the CCVA-measured EC50 values (Fig. 3D, “Model with constrained affinity constants”) while maintaining the affinity constants in a realistic range around their measured values (figs. S14 to S18). We relied on Akaike’s information criterion (see the Supplementary Materials) to select the minimal model that accurately reproduced the CCVA data. This resulted in a model in which extracellular cytokines bind only to the fully assembled receptors IL-7Rα:γc and IL-2Rα:β:γc, and the opposing correlation of the IL-2 and IL-7 EC50 values with IL-2Rα abundance can be readily interpreted as a consequence of the preformation of IL-2R complexes that deprive the IL-7R of the γc chain at high IL-2Rα abundance (Fig. 3E). This prediction of the model could be tested and confirmed experimentally. We sorted effector T cells into two groups, the 35% of cells that had the highest IL-2Rα abundance and the 35% of cells that had the lowest abundance, and analyzed their steady-state association of γc and IL-7Rα by immunoprecipitation and Western blotting analysis (Fig. 3F). This experiment showed that IL-7Rα:γc complexes were present in cells with a low abundance of IL-2Rα, but they were almost totally absent in cells with a high abundance of IL-2Rα.

An extended equilibrium binding model accounts for the effects of IL-2Rα on the EC50 values of other γc family cytokines

IL-7 is only one among the γc family of cytokines whose receptors share some components with IL-2R. Among these, IL-4, IL-7, IL-15, and IL-21 all stimulated a pSTAT5 response in effector T cells, whereas another member, IL-9, did not (fig. S19). With CCVA, we found that the EC50 values of these cells for IL-4, IL-7, IL-15, and IL-21 did not behave identically as a function of IL-2Rα abundance (Fig. 4A). Apart from IL-2, the EC50 of the IL-7 response was the most strongly affected by the abundance of IL-2Rα. For IL-15, we observed a slight dependency, whereas the EC50 values for IL-4 and IL-21 responses were largely unaffected by IL-2Rα abundance. To account for the variation in these behaviors, we extended the simple model selected earlier to include these additional receptors (Fig. 4B), using typical protein abundances defined by flow cytometric calibration (fig. S8).

Fig. 4

Extension of the preassembly model to the full γc family of cytokines. (A to D) Experimental data (A) and model (B to D) of the formation of receptors for IL-2 (blue), IL-7 (red), IL-4 (brown), IL-15 (green), and IL-21 (purple). (A) Correlation of γc-family cytokine EC50 values with IL-2Rα abundance, as measured by CCVA. Error bars represent geometric means ± SD of EC50 values from CCVA of two independent replicate dose responses. Data are representative of four independent experiments. (B) Schematic of the model incorporating all cytokine-specific α chains (see fig. S20 for the different parameters used). The single components are represented by closed circles, whereas complexes are represented by open circles. A given complex sits on the line joining the two reactants involved in the corresponding binding reaction. The signaling cytokine-receptor complexes are colored depending on the cytokine involved. (C) Theoretical dependency of the sensitivity to various cytokines on IL-2Rα abundance for the model shown. (D) Modeled number of preformed cytokine receptors (in the absence of external cytokine) as a function of IL-2Rα abundance. Note that the curve for IL-4 is shifted up slightly to distinguish it from the curve for IL-21.

We modeled the full system with our previously estimated affinity constants related to the formation of IL-2Rs and IL-7Rs, and adjusted the other affinity constants (fig. S20) to reproduce the effects of IL-2Rα abundance on the EC50 values of the responses to these cytokines (Fig. 4C). We found that high-affinity constants for the binding of γc to IL-21R and IL-4Rα were required to recapitulate the independence of the IL-4 and IL-21 EC50 values from IL-2Rα abundance (fig. S20). Consequently, preformation of IL-4Rs and IL-21Rs was largely unaffected by IL-2Rα (Fig. 4D). To our knowledge, affinity constants for the binding of γc to IL-21R and IL-4Rα have not been measured, and our model predicts high values for these affinity constants on cell surfaces. Average EC50 values are tuned through the binding affinity of the preassembled receptor to the cytokine and the abundance of IL-21R or IL-4Rα. Similarly, in the case of IL-15, a high-affinity binding between IL-15Rα and the β:γc dimer (fig. S20) led to a weak correlation between the IL-15 EC50 and IL-2Rα abundance. In our model of the γc cytokine family, the independence of a given cytokine EC50 value from IL-2Rα abundance is achieved when the corresponding receptor α chain has a high affinity for the shared receptor component(s).


Here, we presented a comprehensive methodology to take advantage of the variation of protein abundance within cell populations to quantify the variability in the sensitivity of these cells to different cytokines and to investigate the underlying mechanisms involved. This method relies on the demonstrated ability of flow cytometry to precisely measure the abundance of multiple proteins simultaneously in single cells (Fig. 1). We used this methodology to study receptor signaling by members of the γc chain family of cytokines, and uncovered unexpected relationships between related receptors.

During the antigenic activation of isogenic T cells, the abundance of IL-2Rα varies by up to five orders of magnitude (17). Our costaining experiments verified that these measurements of IL-2Rα precisely defined IL-2Rα abundance in these cells (Fig. 1F and fig. S1). This finding enabled us to study the effect of this broad range of IL-2Rα abundance on cytokine sensitivity, as defined by the EC50 of STAT5 phosphorylation stimulated by the cytokine. In agreement with previous findings, we found that IL-2Rα abundance strongly correlated with IL-2 sensitivity, with an EC50 decreasing from 100 to 1 pM as IL-2Rα abundance increased from 102 to 107 proteins per cell (Fig. 1 and fig. S2). This finding is consistent with the observation by Cantrell and Smith that cells endowed with a higher IL-2Rα abundance proliferate more robustly in vitro than do cells with a lower receptor abundance (22). We observed that this increase in the sensitivity of cells to IL-2 was concomitant with a 20-fold decrease in their sensitivity to IL-7 (Fig. 3). This dependency of IL-7 sensitivity on IL-2Rα abundance could not be explained by an inverse correlation between IL-7Rα and IL-2Rα abundances or by a potential inhibition of IL-7Rα signaling capacity as a function of IL-2Rα abundance within activated effector cells (figs. S3 and S4). Such inverse correlation in receptor abundances has, however, been observed during the transition of T cells from naïve to effector stages and from effector to memory stages in vivo. IL-7Rα is abundant in naïve cells, decreases in effector cells, and is then restored in memory cells and their precursors, whereas IL-2Rα follows the opposite pattern: it is absent in naïve cells, highly abundant in effector cells, and then diminished in memory cells and their precursors (18, 25). Additionally, we found that highly activated effector cells, which had large amounts of IL-2Rα, were inhibited from sensing IL-7, a finding that implies that these cells would rely chiefly on IL-2 for survival. This could reinforce the suppressive effect of IL-2 depletion by T regulatory cells, preventing IL-7 from rescuing effector cells when IL-2 becomes scarce (17, 3739). In this context, the inhibitory effect of IL-2Rα abundance on IL-7 sensitivity represents an additional mechanism to enhance the separation of cytokine specificities across different T cell subtypes, with naïve and memory cells relying chiefly on IL-7, whereas effector cells rely on IL-2.

Addressing the impact of IL-2Rα abundance on IL-7 sensitivity at the computational level was made possible because this system is well characterized at the biophysical level (20, 27, 29, 3136). Binding interactions between receptor subunits and between receptor subunits and cytokines have been quantified (29). However, because many of these measurements have been made in different systems and with varying methodologies—on the surfaces of cells or with soluble fractions—multiple models have been proposed for the functions of these receptors (28). Our CCVA methodology (Fig. 2) enabled us to add detailed quantitative data, across a large range of protein abundance, thereby placing tight constraints on any proposed model of IL-2 and IL-7 signaling. Here, we showed that equilibrium binding models constrained by previous measurements of affinity constants quantitatively reproduced aspects of the dependency of IL-7 and IL-2 sensitivities on IL-2Rα abundance while remaining relatively close to the measured values of the affinity constants (Fig. 3). The switch of sensitivities between IL-2 and IL-7 as a function of IL-2Rα abundance can be readily attributed to the preformation of the full heterotrimeric IL-2R complex, which sequesters the γc chain away from IL-7R. By comparing all possible equilibrium binding models for the formation of signaling cytokine-receptor complexes, we identified a minimal equilibrium binding model with affinity constants constrained around the measured values that accounted quantitatively for the CCVA data. This model was selected with a statistical criterion consisting of a cost for deviation from both biochemically measured parameters and CCVA data, and penalizing for the inclusion of a larger set of parameters. Hence, this selection procedure does not exclude other models with more parameters, representing other paths of receptor formation (35, 40, 41), but more data would be needed to justify their inclusion. The selected model presents the practical advantage of providing a straightforward interpretation of the scaling of IL-7 sensitivity as a function of IL-2Rα abundance. The increased EC50 at high IL-2Rα abundance is a direct consequence of the sequestration of the γc chain in the preformed IL-2R. As a consequence, a cytokine receptor subunit with a high affinity for a shared component (the γc or β chain) should make it difficult for the IL-2R to sequester this component, and the effect of IL-2Rα abundance on the sensitivity to the corresponding cytokine should be weak.

Extending this minimal model to other members of the γc family of cytokines, we recapitulated different measured effects of IL-2Rα abundance on the sensitivities of T cells to cytokines (Fig. 4). Specifically, we found that IL-2Rα abundance reduced the sensitivity of cells to IL-7 and, to a lesser extent, to IL-15, but not to IL-4 or IL-21. In our model, the independence from IL-2Rα abundance of the sensitivity of cells to these two cytokines implies high-affinity binding between the corresponding α chains and the γc chain. The ability of each receptor α chain to tune not only the sensitivity of the receptor to its respective cytokine but also the cross-antagonism of other signals provides a functional separation between γc-using cytokines that matches their respective functions in the context of T cell differentiation (42). IL-7 is dedicated to the homeostasis of naïve and memory cells, but not to the proliferation of effectors, whereas IL-15 affects both effector cells and memory cells (43). Therefore, segregating the homeostatic functions of IL-7 from the proliferative function of IL-2 may be advantageous (24). In contrast, IL-4 and IL-21 are both involved in the differentiation of effector cells into helper subtypes, and cells may be exposed to either or both during differentiation (21). For these two cytokines, cross-antagonism by large amounts of IL-2Rα would thus be functionally counterproductive.

Our CCVA method can be applied to virtually any cellular signaling investigation wherein cellular heterogeneity and signaling can be measured, such as for other cytokines, chemokines, or growth factors (44, 45). The number of signaling regulators that can be simultaneously measured is limited only by the technical constraints of flow cytometry. As the limitations of flow cytometry recede with the introduction of machines capable of recording dozens of parameters simultaneously (10, 11), our method is poised to identify and quantify regulatory effects in signaling networks. This technical framework, accompanied by quantitative modeling efforts, will improve our mechanistic understanding of signal transduction cascades in living cells (2, 3).

Materials and Methods

Reagents and antibodies

Recombinant mouse IL-2, IL-4, IL-7, IL-9, IL-12, IL-21, and TSLP were from eBioscience, and recombinant mouse IL-15 was from R&D Systems. For a list of antibodies used for flow cytometry, cell sorting, immunoprecipitation, and Western blotting, see table S1.

Cell culture media

All in vitro experiments were performed in complemented RPMI medium (prepared by the media facility at Memorial Sloan-Kettering Cancer Center), which consists of RPMI 1640 supplemented with 10% heat-inactivated fetal bovine serum, 2 mM l-glutamine, 10 mM Hepes (pH 7.4), 0.1 mM nonessential amino acids, 1 mM sodium pyruvate, penicillin (100 μg/ml), streptomycin (100 μg/ml), and 50 μM β-mercaptoethanol. All cell cultures were maintained in an incubator at 37°C with 5% CO2.

Cell culture

All cells used in this work were 5C.C7 TCR–transgenic, Rag2−/− IL-2−/− CD4+ T cells, isolated and cultured as described previously (13). Briefly, cells were incubated with irradiated B10A CD3−/− splenocytes pulsed with 1 μM MCC K5 peptide: ANERADLIAYFKAATKF (GenScript). IL-2 (10 nM) was given at 1 and 4 days after peptide stimulation, and then every 3 days for the maintenance of cultures. Cells were used between day 5 and day 20 of culture, unless otherwise noted. Live cells were purified before analysis by separation on a Ficoll gradient (GE).

Cytokine titrations

Before measuring cytokine responses, bound IL-2 was stripped from the cell surface by a 2-min incubation in 0.1 M glycine buffer equilibrated at pH 4.0. After two washes in RPMI, cells were rested for 30 min at 37°C (46). Cytokine titrations were prepared in 96-well V-bottom plates on ice. Cells (105 per well) were added to cytokines and incubated for 10 min at 37°C before being processed for flow cytometric analysis.

Flow cytometry

Cells were prepared for single-cell staining by following protocols optimized by the Nolan group (47). Briefly, cells were fixed for 10 min on ice in 1.6% paraformaldehyde followed by permeabilization on ice in 90% methanol. Cells were then washed twice in fluorescence-activated cell sorting (FACS) buffer [phosphate-buffered saline (PBS) with 4% fetal calf serum and 0.1% sodium azide] and stained intracellularly with antibody specific for pSTAT5, followed by incubation with labeled secondary antibody and conjugated antibodies specific for receptors (CD4 and IL-2Rα). All staining was performed at room temperature in FACS buffer. Before flow cytometric acquisition, cells were stained with 5 μM 4′,6-diamidino-2-phenylindole (DAPI; Sigma-Aldrich) to exclude dividing (G2) and apoptotic (sub-G1) cells. Flow cytometry acquisition was performed on an LSR-II flow cytometer (BD Biosciences).

Receptor abundance calibration

A standard curve correlating fluorescence as measured by flow cytometry to the number of phycoerythrin (PE) molecules was calculated with Quantum R-PE molecules of equivalent soluble fluorochrome (MESF) calibration beads (Bangs Laboratories). “Corrected median fluorescence intensity” (cMedFI) values were calculated by subtracting the median fluorescence of unstained beads from the median fluorescence value of each PE standard. A linear regression of log cMedFI to log PE bead MESF values was used to establish a standard curve. To calibrate typical abundances of receptor chains, we stained live cells with saturating concentrations (as established by titration) of the respective PE-labeled antibodies or a nonspecific PE-labeled rat immunoglobulin G (IgG) control. Corrected median fluorescence was calculated by subtracting the median fluorescence of the isotype control from the median fluorescence of each antibody staining. The calibration curve was used directly to calculate the numbers of molecules from the corrected median fluorescence (fig. S8) using a ratio of one PE molecule per antibody, as specified by the manufacturer (eBioscience or BD Biosciences). For calibration of IL-2Rα abundance in ScatterSlice data, the calibration curve was applied to the fluorescence value of each bin.

Population sensitivity fitting

We determined the EC50 values for given cellular populations from flow cytometry data by fitting the dose response for the geometric mean of pSTAT5 as a function of the cytokine concentration with the following equation:pSTAT5([cytokine])=pSTAT5low+pSTAT5highpSTAT5low1+EC50[cytokine],where pSTAT5low and pSTAT5high are the baseline and the plateau of pSTAT5 signals. Fits were performed with Prism (GraphPad Software). The fitting algorithm uses the Hessian matrix (48) to estimate errors on the fitted parameters.

Fluorescence-activated cell sorting

Cells were first purified by Ficoll and then incubated with PE-conjugated anti-CD25 antibody (clone PC61.5, Miltenyi Biotec) and allophycocyanin-conjugated anti-CD4 antibody on ice for 30 min. Before sorting, 5 μM DAPI was added to exclude dead cells. Sorting was performed on a FACSAria flow cytometer (BD Biosciences). Sorted cells were kept on ice until further processing.

Western blotting analysis

Equal numbers of sorted cells were rinsed in PBS and then lysed in buffer containing 25 mM tris-HCl, 1 mM EDTA, 130 mM NaCl, and 1% NP-40, supplemented with protease inhibitor cocktail (Sigma-Aldrich), which resulted in final concentrations of 1.04 mM AEBSF [4-(2-aminoethyl)-benzenesulfonyl fluoride hydrochloride], 800 μM aprotinin, 40 μM bestatin, 14 μM E-64, 20 μM leupeptin, and 15 μM pepstatin A. Lysates were mixed with NuPAGE LDS loading buffer with 2.5% β-mercaptoethanol and were incubated at 95°C for 2 min before being loaded and resolved on NuPAGE Novex 4 to 12% bis-tris gels (all products from Life Technologies). Gels were transferred onto Immobilon membranes (Millipore) before being incubated with primary antibody against IL-7Rα followed by horseradish peroxidase–conjugated anti-rabbit IgG secondary antibody in tris-buffered saline with 0.05% Tween 20 and 2% bovine serum albumin.

Immunoprecipitation of the γc chain from sorted cells

Cultured T cells (5 × 108) were stripped of cytokines in 5 ml of 0.1 M glycine buffer equilibrated at pH 4.0 to remove any remaining IL-2 from the culture. After two washes in RPMI, cells were rested for 30 min at 37°C. Cells were centrifuged and then resuspended in 10 ml of fresh RPMI for staining and were processed for FACS analysis as described earlier. The 35% of cells that had either the greatest or the least abundance of IL-2Rα were collected. After cell sorting, immunoprecipitations were performed with the Dynabeads coimmunoprecipitation kit (Life Technologies) following the manufacturer’s protocol, with the following modifications. Cells were centrifuged, weighed, and suspended in a 13:1 ratio (μl:mg) of Dynabeads lysis buffer [“buffering salts” (pH 7.4), 110 mM potassium acetate, 0.5% Triton X-100, 100 μM sodium chloride]. After incubation on ice for 20 min, cell lysates were passed through a 27-gauge syringe 10 times on ice. Immunoprecipitations were performed according to the manufacturer’s instructions, with Dynabeads freshly conjugated to anti-γc antibody (clone H-300). Western blotting was then performed as described earlier.

ScatterSlice analysis

Initial analysis of flow cytometry data was performed with FlowJo software (TreeStar). Data corresponding to single CD4+ cells with 2N DNA content were identified and exported as text files for analysis by our processing R program ScatterSlice (49). The software is made available as supplementary material for this paper. Details of the algorithm are as follows. Text files are imported into R and then divided into user-specified bins. Within each bin, a three-parameter Hill equation, which fits the base, amplitude, and EC50, is fit by minimizing the sum across all concentrations in a dose response (from i = 1 to Ndoses) of the squared difference between the log of each dose’s geometric mean fluorescence intensity of the response channel (pSTAT5 in our studies) (ziobs) minus the log of the three-parameter Hill equation (fitting base, amplitude, and EC50), normalized by multiplying by the number of cells in that dose (Ni) divided by the phospho-response variance (σi2). Errors on fitted parameters are estimated from the Hessian matrix (48). Tables of fitting results (values for base, amplitude, EC50, and their errors for different amounts of receptors) are presented as color maps or line graphs and can be further analyzed with R or exported as text files for analysis with other software.

Supplementary Materials

Computing cytokine sensitivity from affinity constants

Maximum a posteriori estimation of the affinity constants

Model enumeration and selection

Error estimation

Fig. S1. Flow cytometry provides accurate measures of relative protein abundance.

Fig. S2. CCVA analysis of cell populations from Fig. 1.

Fig. S3. IL-7Rα abundance does not correlate with variation in IL-7 EC50.

Fig. S4. Responsiveness to TSLP demonstrates the signaling functionality of IL-7Rα.

Fig. S5. Kinetics of the chemical equilibrium for the pSTAT5 response.

Fig. S6. The pSTAT5 response is linear with the number of bound IL-7 molecules per cell.

Fig. S7. Cytokine depletion is negligible and does not affect the IL-2 or IL-7 EC50 values.

Fig. S8. Methodology for absolute receptor quantification.

Fig. S9. Estimation of the set of “measured” affinity constants, Kexp.

Fig. S10. A nondirect influence of α2 on γ cannot reproduce the dependency of IL-7 sensitivity on α2.

Fig. S11. Indexing of all relevant models.

Fig. S12. The experimentally determined affinity constants cannot quantitatively account for the effect of IL-2Rα on IL-2 and IL-7 EC50 values.

Fig. S13. Fitting without constraints on the affinity constants leads to great divergence from their measured values.

Fig. S14. Determination of the optimal value of ρ as a parameter constraint for model fitting.

Fig. S15. Model selection with corrected AIC.

Fig. S16. Determination of the 95% uncertainty intervals for the parameters from the Profile Likelihood (PL).

Fig. S17. Two-dimensional sections of the posterior distribution for the equilibrium constants for ρ → ∞ (unconstrained fit).

Fig. S18. Two-dimensional sections of the posterior distribution for the equilibrium constants for ρ = 1 (constrained fit).

Fig. S19. pSTAT5 dose responses for γc family members show measurable responses and EC50 values.

Fig. S20. Set of parameters for the extended γc family model.

Table S1. List of antibodies used in this study.

ScatterSlice software and instruction manual

References and Notes

Acknowledgments: We thank J. Sethna, K. C. Garcia, K. Tkach, N. Altan-Bonnet, M. Vergassola, and T. Rose for useful comments and discussion. Funding: This work was supported by NSF grant 0848030 and NIH grants R01-AI083408, U54-CA143798, and T32-AIO7621 (J.W.C.). Author contributions: J.W.C., G.V., and G.A.-B. designed the research; J.W.C. and V.K. performed the experiments; G.V. performed the computational modeling; J.W.C. and O.E.D. designed and coded ScatterSlice; and J.W.C., G.V., and G.A.-B. wrote the manuscript. Competing interests: The authors declare that they have no competing interests.

Stay Connected to Science Signaling

Navigate This Article