Quantitative Network Signal Combinations Downstream of TCR Activation Can Predict IL-2 Production Response¹

Signaling through the TCR leading to cell decision processes is determined in part by the affinity of the TCR for its ligand. Affinities of TCR/ligand interactions can differ by at least three orders of magnitude (1) and can correspondingly elicit vastly different cellular responses. In mature T cells, high-affinity TCR-peptide MHC (pMHC)⁴ interactions can lead to T cell activation, whereas low-affinity interactions are important for cellular homeostasis and may be involved in antigenic signal amplification (2, 3, 4). Activation of T cells by altered peptide ligands (APL) can also induce partial activation, hyperstimulation, and anergy (5). Elucidating the complex molecular details underlying ligand affinity discrimination by T cells is central to our understanding of both T cell-mediated immunity to pathogens and immune dysfunctions, such as autoimmunity and aberrant immunosurveillance.

Cell fate responses to TCR ligation can vary with single residue changes occurring on the peptide at the ligand-receptor interface. Upon TCR engagement, changes in the ITAM phosphorylation patterns (6, 7) activate multiple interconnected signaling pathways (8, 9) (illustrated in Fig. 1 A). Propagation of stimuli from the receptor results in T cell activation and nuclear localization of AP-1, NFAT, NF-κB, and Oct-1 transcription factors. One early marker of T cell activation is IL-2 production. AP-1, NFAT, NF-κB, and Oct-1 transcription factors are all known to bind to the IL-2 promoter region (10).

A growing body of work has focused on proximal TCR signaling events to individual downstream pathways to elucidate mechanisms of ligand affinity discrimination (3). However, there exists little understanding about how the dynamic state of the multipathway network activated by TCR signaling, as well as by costimulatory cues, may explicitly encode the information crucial for the ultimate cell behavioral responses governed discriminatorily by TCR ligation events. It is worth noting that recent studies in other areas of cell biology have ascertained that quantitative combinations of nodes in multiple downstream pathways across a signaling network distal to initiating receptor activation integratively represent the receptor-mediated information governing cell phenotypic behavior (11, 12, 13). In this study, we propose that in the context of TCR stimulation, activation states quantitatively encoded by the signaling network define regulatory aspects for a T cell cytokine production response.

To test this hypothesis, a T cell hybridoma line was stimulated with APL and a quantitative data set of both signaling dynamics and cellular response data was generated across the full set of APL treatments. IL-2 production was chosen as a relevant cell functional response, and signaling network nodes were defined on the basis of prior literature knowledge. Because our integrative systems approach of quantitative measurement of multiple protein signaling components across multiple cell treatment cue conditions requires relatively large cell numbers, this initial proof-of-concept study was conducted with a T cell hybridoma cell line. To gain predictive capability and biological insights from this multivariate data compendium, we used a high-level computational modeling approach (14, 15), partial least squares regression (PLSR) analysis, to understand the relationship between ligand avidity, downstream signaling dynamics, and IL-2 production. PLSR analysis is a multivariate regression technique that can relate a data set of predictor variables (signaling measurements) to a response variable (IL-2 production) without prior knowledge of—or assumptions about—network structure. This method of linear multivariate statistical modeling has been used successfully to understand signal-response relationships in complex cell decision processes such as apoptosis (12, 16, 17). In this study, we demonstrate through a model hybridoma system that this approach can be used successfully to understand the relationship between the strength of ligand-receptor interaction and consequent signaling pathways triggered by TCR engagement.

The 1B6 T cell hybridoma was generated in the laboratory of Vijay Kuchroo by fusion of the 1B6 T cell clone (18) with the TCR-negative fusion partner BW1100. Following fusion, positive hybrids were identified by measuring growth inhibition in the presence APCs and L144 peptide vs control peptide. Positive cells were subcloned using limiting dilution at least twice.

DAS and 1B6 hybridomas were cultured in RPMI 1640 with glutamine (VWR), 10% FBS (HyClone), penicillin (100 U/ml), streptomycin (100 μg/ml), nonessential amino acids (1×), HEPES (10 mM), and sodium pyruvate (1 mM) (all VWR).

DAS cells (B7.1-positive fibroblasts transfected with I-A^s) (19) were grown to confluence on Nunc-coated omni trays (Nunclon). Q144 (HSLGKQLGHPDKF), Y144 (HSLGKYLGHPDKF), L144 (HSLGKLLGHPDKF) peptides (Anaspec) were added at concentrations indicated to DAS medium 4–12 h before addition of T cells; 10⁷ 1B6 hybridoma cells were added to DAS plates, which were centrifuged at 1,000 rpm for 1 min, and incubated at 37°C. Supernatant was transferred to 4 vol of ice-cold PBS and centrifuged at 3,000 rpm for 1 min. The cells were lysed in buffer (20) modified to 2% Triton X-100 on ice for 20 min. After centrifuging at 14,000 rpm for 10 min at 4°C, the soluble fraction of lysate was collected (cytosolic extract), and protein concentration was determined using a MicroBCA protein assay kit (Pierce). The pellet was resuspended in 420 mM NaCl, 20 mM HEPES (pH 7.9), 1 mM EGTA (pH 8.0), 1 mM EDTA (pH 8.0), 1 mM sodium orthovanadate, 1 mM DTT, 1 mM PMSF, 10 μg/ml aprotinin, 10 μg/ml leupeptin, 1 μg/ml pepstatin, and 1 μg/ml microcystin LR. The pellet was sonicated for 1 min followed by 3 min on ice, repeated 6 times, followed by centrifugation at 14,000 rpm for 30 min at 4°C.

PD098059, LY294002 (DMSO stock; Calbiochem), cyclohexamide (EtOH stock; Sigma-Aldrich) was diluted to working concentration in medium. 1B6 cells were preincubated in PD98059 and/or LY294002, and cyclohexamide for 30 min before stimulation with peptide-loaded DAS cells.

1B6 cells were stimulated for 0, 1, 2, 3, or 4 h with 4 μg/ml L144 peptide-loaded DAS cells. Total RNA was isolated and purified using TRIzol reagent (Invitrogen Life Technologies) and RNeasy kit (Qiagen) according to manufacturer’s instructions. cDNA was generated from 400 ng of DNase-treated total RNA using random hexamer primers and Ominscript reverse transcriptase (Qiagen). Real-time PCR was performed using a CHROMO4 detector (MJ Research) on 0.5 μl cDNA in 25 μl of reaction volume containing SYBR Green universal PCR master mix (Qiagen) using primers for IL-2 and β-actin (Operon Biotechnologies) and thermal cycling conditions according to Aube et al. (21). To extrapolate Ct values, a manual threshold was set in the Opticon Monitor 3 software package (MJ Research) and kept consistent across all samples. Following real-time PCR, a melting curve was performed from 50–90°C to insure the uniqueness for product and absence of primer dimers. All samples are normalized to β-actin mRNA levels within sample.

For SDS-PAGE, lysates were diluted to either 2 mg/ml or 4 mg/ml, reducing Laemmli-SDS sample buffer was added, and samples were boiled for 10 min. Forty or 80 μg of lysate was loaded per well. Each sample was loaded in triplicate wells and for normalization purposes, each gel contained duplicate lanes of a positive control sample. After electrophoresis, samples were transferred to polyvinylidene difluoride membrane and blocked with 3% BSA in TBST. Blots were visualized and quantified using a Kodak Image Station 1000. The net intensity of individual bands was determined and normalized to the net intensity of the average of the positive control bands. Triplicates were then averaged and fold activation determined by normalization to the zero time point. Anti-phospho-p44/42 (Thr²⁰²/Tyr²⁰⁴), anti-phospho-SAPK/JNK (Thr¹⁸³/Tyr¹⁸⁵), anti-phospho-p38 (Thr¹⁸⁰/Tyr¹⁸²), anti-phospho-Akt (Ser⁴⁷³) were obtained from Cell Signaling Technologies, and anti-rabbit HRP was purchased from Santa Cruz Biotechnology.

MK2 and IκB kinase (IKK)αβ activities were determined by a custom multiplex kinase assay (20) with the following modifications: For the MK2 assay, protein A microtiter strips (Pierce) were coated with anti-MK2 (StressGen Biotechnologies), and the reaction was allowed to proceed for 1 h. For the IKKαβ assay, 500 μg of total lysate protein was incubated in the Ab-coated wells overnight. Data were normalized by the value quantified for the zero time point for fold activity.

NFAT phosphorylation and localization was measured by resolving 22 μg of sample protein (boiled for 10 min in sample buffer) on 5% polyacrylamide gel, alternating cytosolic and nuclear fractions in the lanes. Membranes were blocked with 3% BSA in TBST followed by incubation with NFATc2 Ab (Santa Cruz Biotechnology) at a 1/500 dilution. Ratios of dephosphorylated bands to phosphorylated bands in each lane were determined.

Supernatant was collected after 4 h and analyzed using quantitative sandwich ELISA for mouse IL-2 (BD Biosciences).

The compiled data matrix was transformed by the equation y = log10(x + 1) where x is a measured value, to account for a skewed distribution of experimentally measured values. The column transformed data were then mean-centered and scaled by unit variance. PLSR analysis was performed in SIMCA-P (Umetrics). The data set was divided into two matrices: Y ε R^15x1 consisting of IL-2 measurements, and X ε R^15x57 denoting the measured protein signals. The PLSR algorithm successively finds linear combinations of the predictor variables that have maximum covariance with the response variables, so that each linear function is uncorrelated with previous ones (22). PLSR can accommodate data sets that are not fully complete (i.e., yielding matrices that are not of full rank), provided the missing values are randomly distributed. For more detailed description on PLSR modeling we refer to a previous explication (22).

In T cells containing the 1B6 TCR, the three altered peptides we use have been shown to affect proliferation, cytokine production, and anergy hierarchically (18, 23, 24). The APLs are based on the autoantigen myelin proteolipid protein peptide 139–151 with a single amino acid substitution at position 144 from tryptophan to glutamine, tyrosine or leucine. All three peptides bind MHC class II with equivalent affinities (18). We assayed apoptosis induced by different concentrations and avidities of APLs in T cell hybridomas (25) that express the 1B6 TCR. DAS APCs were pulsed with proteolipid protein APLs Q144, Y144, and L144 for 4–12 h before addition of 1B6 cells. After 24 h, apoptosis of 1B6 cells was measured via propidium iodide incorporation (Fig. 1,B). Prior experiments with annexin V and propidium iodide double staining confirmed that cell death through peptide stimulation occurs by apoptosis (data not shown). 1B6 cells stimulated with L144 peptide showed >85% cell death at both peptide concentrations measured, the Y144 peptide induced 53% cell death at 4 μg/ml, whereas the Q144 peptide at 4 μg/ml induced cell death at levels similar to the no peptide control (Fig. 1 B).

IL-2 secretion has been previously reported in the literature to serve as an appropriate marker of T cell activation in T cell hybridoma lines; this response is specific to TCR stimulation, with quantitative differences in IL-2 levels observed with different peptide binding (26, 27). In the 1B6 hybridoma cells, the levels of IL-2 production from stimulation with the various peptides correlated with the apoptotic responses (Fig. 1 C). IL-2 production was measured after 4 h of peptide stimulation. High levels of IL-2 were present in the medium of 1B6 cells stimulated with L144 peptide and low levels of IL-2 were present in the medium of cells stimulated with Y144. IL-2 in the supernatant was below the detection limit in cells stimulated with Q144 peptide.

To determine whether IL-2 production is transcriptionally regulated following peptide stimulation, IL-2 mRNA levels were measured using quantitative RT-PCR (Table I). 1B6 hybridomas stimulated with L144 peptide showed increased levels of IL-2 mRNA after 1 h of stimulation and remained elevated for the duration of the experiment. To determine whether IL-2 production is regulated at the level of protein translation following peptide stimulation, 1B6 cells were stimulated with L144 peptide in the presence of cyclohexamide (Table II). IL-2 production in cyclohexamide-treated samples is reduced to levels comparable to no-stimulation controls. These results are consistent with IL-2 production in 1B6 hybridomas being regulated at the level of both transcription and translation following peptide stimulation.

To learn how the avidity of the TCR for the pMHC might be encoded within these signaling pathways, we measured the activation levels of key signaling molecules known to be important for IL-2 production in T cells (highlighted molecules in Fig. 1,A). Of the IL-2 promoter elements, we measured the nuclear translocation of one directly (NFAT), and measured activation/phosphorylation of kinases upstream of others (IKKαβ for NF-κB, JNK, ERK1/2, and p38/MK2 as activators of AP-1). To generate a self-consistent, multivariate data set (17, 20), all protein activation measurements (Fig. 1,A) were obtained from a single lysate containing in excess of 40 × 10⁶ cells per stimulation condition and time point. Lysates generated from the APCs after coculture with the 1B6 T cells show no change in the activation state of signaling molecules measured (data not shown), confirming that the dynamic time courses represent changes in T cell signaling states. The end point cellular response measurement was the level of IL-2 production in the supernatant collected from the plate at 4 h. In total, we measured five peptide stimulation conditions and one no peptide control condition for 11 predictor signals at six time points, and one response variable (IL-2) at the final 4-h time point. To account for biological variation within the data set, each stimulation condition was replicated on three different days and triplicate experiments were entered into the model as separate observations (Fig. 2,A). In a few cases, we were unable to accurately quantify signals on a lysate (Fig. 2 A, gray squares). These data points were not excluded for statistical reasons, rather the bands on the membrane were unquantifiable for technical reasons. The final data matrix contained ∼1,200 independent measurements. To our knowledge, this is the first time that such a broadly distributed array of signals across multiple pathways has been analyzed systematically and concurrently for a clonal TCR line activated by peptide ligand stimulation.

Two illustrative examples of the dynamic time courses of these signaling measurements are shown in Fig. 2,B. The replicate time course measurements for ERK2 and Akt are each compared across the three different peptides all at a concentration of 4 μg/ml. The curves in these plots follow similar dynamics for the various peptide conditions, differentiated mainly by quantitative magnitudes of the various peaks. Perusal of the other signaling nodes in Fig. 2 A indicates that the left-to-right time courses for a given node seem to be at least roughly similar in profile going down the corresponding peptide treatment column. Thus, within the context of this experimental system, the effects of the various APL appear to be mainly on the relative magnitudes of the dynamic activation profiles of the signaling nodes sampled here. This finding may reflect in explicit, quantitative fashion the fact that the signaling measurements examined here largely represent nodes downstream of the TCR triggering processes. Qualitatively different signatures of membrane-proximal signaling markers have been observed upon T cell stimulation by APLs; e.g., differences in patterns of TCR ζ-chain phosphorylation have been found (28). Such differences could be expected to be less prominent for downstream signaling nodes if the APLs considered all stimulate membrane-proximal signaling that exceeds a threshold.

Our objective in this work was to explore the hypothesis that the strength of the TCR stimulus is encoded within the signaling dynamics of the network downstream of TCR activation. If the levels of IL-2 production are indeed encoded within the signaling dynamics then it is feasible to construct a model capable of predicting the IL-2 production data in Fig. 2 A from quantitative combinations of the signaling data. PLSR analysis is particularly useful for the analysis of biological data because it can handle noisy, collinear, and incomplete variables (16). This technique permits determination of the dependence of response variables (IL-2 production) on quantitative combinations of a key subset of predictor variables (the dynamic signals). These combinations are termed principal components and can be thought of as “weighted” signals, whereas the predicted response is made in terms of “loadings” applied to the principal components.

A PLSR model built on the above described signaling measurements and response variables converged at two principal components. The two-principal component model is capable of fitting IL-2 concentrations from the input data with a line y = 1.19x + 0.06 and a regression coefficient (defined as R²Y; see Materials and Methods) of 0.97. As mentioned, these principal components are linear combinations of predictor variables, posing a potential challenge for highly connected biological pathways where nonlinearities may be present in the system. We examined the potential benefit of including nonlinear terms to the model, as described previously (16), but found that these additional terms did not enhance the model’s ability to explain data variance; in fact, they decreased the regression coefficient by amplifying noise in the data matrix (data not shown). This suggests that a linear PLSR model can describe the highly nonlinear TCR pathway, from the activation state of signaling molecules to an IL-2 output. A linear PLSR model has previously been shown to describe the apoptotic response network (12, 16, 17).

To determine how individual samples (defined as all measurements for a replicate of a given stimulation condition) influence the model’s IL-2 prediction, samples from different experimental conditions can be mapped onto the two principal components via their scores. Individual samples that fall close to the IL-2 response are more heavily weighted by the model for IL-2 prediction and samples that cluster together have similar information content. When individual samples of the TCR response data set are mapped onto the two principal components via their scores, IL-2 producers can be discriminated from non-IL-2 producers on the first principal component (Fig. 3,A). L144 conditions, the substantive IL-2 producers, scored almost exclusively positive on the first principal component while non-IL-2 producers scored negative on the first principal component; this discrimination was found to be highly significant (p < 0.001) by a two-tailed Welch’s t test (df = 13). In contrast, the second principal component loadings are not statistically different between these two groups, so we can infer that the quantitative combination of signals captured within the first principal component represents an integrative “on-off switch” predominantly governing IL-2 production. At the same time, the second principal component was found necessary for generating an optimal model, suggesting that additional information does reside within this component—perhaps “fine-tuning” the magnitude of IL-2 production. Thus, although both principal components one and two work together to specify IL-2 production levels, Fig. 3 A shows a nontrivial distinction between their respective contributions.

In a manner analogous to determining the influence of an individual sample on the PLSR model, one can determine the influence of individual predictor signaling variables on IL-2 prediction. For instance, “Akt at 30 min” defines one such variable. Signals that map close to the IL-2 response position are more highly correlated with IL-2 than are signals that are close to the origin. Signals with large values on either axis are heavily weighted in the respective principal component. To determine how individual predictor signaling variables influence the IL-2 prediction, these variable weights were plotted onto the two principal components based on their contributions to the principal component (Fig. 3, B and C). Signals with large weighted contributions influenced the model to a greater extent than signals that were close to the origin, whereas signals with negative values were inversely related to the response variable (Fig. 3, B and C; the complete mapping can be obtained by contacting the corresponding author). Both ERK1 and ERK2 were heavily weighted at multiple time points (Fig. 3,B). This is consistent with expectations, because ERK is thought to be a major contributor to cellular outcome in response to TCR signaling. More surprisingly, Akt also emerged as even more strongly positively weighted on both the first and the second principal components (Fig. 3 C). Although Akt has been shown to be activated downstream of the TCR and is known to contribute to costimulation, a precise functional role for Akt in TCR signaling has not been reported. Because of the unexpected influence of Akt on both components in the model, we asked whether costimulation significantly contributed to activation of the sampled signaling molecules. We found that the 1B6 cell line does not exhibit noticeable changes in IL-2 when potential CD28/CTLA-4 interactions are blocked with a CTLA-4/Fc fusion protein, nor are changes in anti-CD3-induced apoptosis observed in the presence of anti-CD28 (data not shown), suggesting that the observed influence of Akt in the model is due to primary TCR signaling.

Finally, the weighting of JNK was found to be distributed widely and relatively mildly across both principal components; moreover, its relationship to IL-2 production appeared to roughly colocalize with that of the dephosphorylation-to-phosphorylation ratio of cytosolic and nuclear NFAT (data not shown). This colocalization may reflect an established molecular mechanism that JNK can directly activate NFAT-driven transcription (29), and the widely distributed weighting may correspond to the complex modulation of IL-2 production by the dynamic shuttling of NFAT between cytosolic and nuclear compartments.

We sought to determine whether the biological information captured in the model is capable of predicting IL-2 production for a new stimulation condition. 1B6 T cells were treated with a new condition, 0.04 μg/ml L144 peptide, and the protein signaling variables were input into the model (Fig. 4,A). The model was then used to predict IL-2 production based upon the signal/response relationships derived from the trained model. Fig. 4,B shows the linear regression fit between experimentally measured IL-2 values and the fitted model IL-2 values for the trained model described by Fig. 3,A (black, blue, green, and red symbols). On the same graph, the experimentally measured IL-2 values vs the predicted IL-2 values (yellow symbols) are plotted for the data in Fig. 4 A. The model performs very successfully for two of the samples, predicting IL-2 concentrations within 2% and 18% of the experimentally measured values. For the third sample, model performance was less tight but nonetheless predicted IL-2 concentration to within a factor of 2.

To investigate whether a model trained on measurements of a single signaling protein (with full time course; Fig. 5, A and B) or at a single time point (across the full set of proteins; Fig. 5, C and D) could adequately predict IL-2 production, we constructed different models from restricted subsets of our data. We evaluated the quality of each model by three metrics: the fit to IL-2 (R²Y; a value near 1 is desired), the fraction of variation in IL-2 measured by cross validation (Q²Y; a value near 1 is also desired), and the sum of errors in the three a priori predictions (ΣDmodY; a value near zero is desired).

The optimal full model exhibited a value of 0.94 for R²Y, 0.57 for Q²Y, and 2.04 for ΣDmodY. By comparison, the best-performing single-component model arose from that constructed on ERK2, with slightly improved Q²Y and ΣDModY of 0.6 and 1.58, respectively, but a far less optimal R²Y of 0.81 (Fig. 5, A and B). All other single-protein models performed worse than the ERK2 model. Although Akt had been found to be a primary informer in the full model, an Akt-only model performed poorly by itself, producing an R²Y of 0.62. We can conclude then that a satisfactory model relating signals downstream of TCR activation to IL-2 production requires a quantitative combination of measurements across multiple pathways. This conclusion is reasonable in light of the necessity for multiple transcription factors regulated by the diverse pathways to bind in combinatorial fashion to the IL-2 promoter region.

To address the inclusion of time course data, models were also constructed across the signaling proteins at a single time point (Fig. 5, C and D). The large number of t = 30-min variables identified to be heavily weighted in the full model suggested that a model using only 30-min measurements of the 11 signaling proteins might perform well. An R²Y = 0.90 was indeed found for a model containing only 30-min predictor variables, but the predictive ability of the model was low as evidenced by the large value of ΣDmodY = 7.5. Models built using other single time points performed better than the t = 30-min model for ΣDmodY but at the expense of a lower R²Y. As with the single protein models, inclusion of measurements at multiple time points was required for optimal model prediction. These simulation results indicate that predictive information encoded in the various signaling pathways may differ temporally in kinetic rates of propagation of information relevant to IL-2 production. Therefore, we conclude that to most effectively capture the critical features of network operation, dynamic measurements and measurements across multiple signaling pathways are required.

We have so far demonstrated that a univariate model (i.e., one in which only a single signal, such as ERK or Akt, is included by itself) cannot effectively predict the IL-2 response. But, because most typically therapeutic agents are aimed at individual signaling components, we can ask whether our multivariate model (in which the full complement of multiple pathway information is included) can successfully predict the net effect of a single-target (e.g., ERK or Akt) inhibitory drug on the ultimate cell response. Because ERK has a well-established role in propagating TCR signals and is heavily weighted on the first principal component, we examined the predicted effect of inhibiting ERK on IL-2 production. In addition, we examined the predicted effect of inhibiting Akt on IL-2 production because the optimal model relied especially heavily on Akt in both the first and the second principal components. Simulated time courses were generated by computationally reducing the levels of either ERK or Akt phosphorylation by 25, 50, or 75% across all time points while keeping all other predictor variables unchanged. We modified four signaling time courses: two L144 4 μg/ml observations used in the training set and two L144 0.04 μg/ml observations used in the a priori prediction. The PLSR model was then used to predict IL-2 production for the new observation based upon the altered ERK or Akt signals (Fig. 6, A and B). Reduced IL-2 production levels were predicted to be tightly correlated with both reduced ERK and reduced Akt phosphorylation. The predictions indicate that Akt inhibition will affect IL-2 production levels to a greater extent than ERK inhibition. The effect on IL-2 production of ERK inhibition can be described by a slope of −0.28, and Akt inhibition can be described by a slope of −0.83. Furthermore, the model predicted an unanticipated negative synergy between the ERK and Akt contributions to governing IL-2 production: reducing ERK and Akt levels by the same degrees simultaneously is predicted to yield a smaller effect on IL-2 than the sum of the individual ERK and Akt reductions (Fig. 6 C).

To test these a priori model predictions directly, we experimentally measured IL-2 levels under conditions of ERK and Akt inhibition using small molecule inhibitors and compared these values to the corresponding PLSR model predictions. An experimental dose-response curve was created for inhibition of MEK/ERK by PD98059 and PI3K/Akt by LY294002 (data not shown). 1B6 cells were preincubated with either PD98059, LY294002, or a combination of the two, followed by a 30-min stimulation with 4 μg/ml L144 peptide. Cell lysates were analyzed for phosphorylation state of the inhibition target. PD98059 was found to inhibit ERK phosphorylation with an IC₅₀ of 0.7 μM for ERK1 and 3.1 μM for ERK2. Akt phosphorylation was inhibited by LY294002 with an IC₅₀ of 1.3 μM. To assess the effect of reduced ERK activity on IL-2 production, we inhibited ERK phosphorylation levels by concentrations 0.6, 1.7, and 5 μM based upon the above determined dose response and measured 4 h IL-2 production under these conditions (Fig. 6, A and B). As anticipated, we observed a tight correlation between measured IL-2 levels and ERK inhibition levels, described by a slope of −0.58. Similar experiments using the Akt inhibitor LY294002 at 25, 50, and 75% inhibition (0.5, 1.3, and 5.2 μM, respectively) validated the model-predicted dependence of IL-2 production on Akt phosphorylation levels, with an experimental slope of −1.22. Moreover, as predicted by the model, the slope of Akt inhibition is steeper than that of ERK inhibition. Also as predicted, the decrease in IL-2 due to the combination of ERK and Akt inhibitors was less than the sum of individual inhibitions (Fig. 6,C). Finally, we note that analogous but univariate models based on only an individual signal such as ERK or Akt, as statistically characterized in Fig. 5, are not as effective as the full multivariate model in predicting these drug inhibition effects on IL-2 (see Fig. 7).

By examining the relationship between the activation state of the multipathway signaling network downstream of TCR-pMHC engagement and associated levels of IL-2 production, for APL possessing diverse TCR binding avidities, in an integrative systems manner, we are able here to provide new insights concerning network-wide information processing from TCR/ligand binding to ultimate cell functional responses. Previous work on TCR ligand discrimination has focused on molecules proximal to the TCR; by moving distal from the receptor we have gained a better understanding of how the consequent dynamic state of the TCR signaling network integrates molecular activation dynamics to govern an important T cell phenotypic response. Capturing the multiple activation states generated by a signaling network through an appropriate computational model has been recently shown to be an especially effective approach to understanding complex biological information processing in other cell biology contexts (11, 12, 13).

We chose a TCR/APL system that is well suited to probe the relationship between peptide ligand strength and signaling dynamics due to a multiple-ligand repertoire that can alter the magnitude of an IL-2 response. Another interesting question would have been to compare the signaling difference between agonist and antagonist peptides. However, with the 1B6 system we were limited to studying agonist peptides because antagonist peptides for this system have not been identified. The cytokine output was selected because of the well-defined IL-2 promoter region and the detailed knowledge of its transcription factor recruitment. An attractive next step would be to determine how primary cells integrate signals from costimulatory molecules to regulate IL-2 production and other facets of T cell activation. A limitation of the current study is the use of a hybridoma system capable of an IL-2 production response in the absence of costimulatory signals, but we believe that our successful proof-of-concept study strongly motivates a next attempt to move our approach into primary cell studies. Another beckoning line of inquiry is understanding how the more well-understood TCR-proximal events propagate to quantitatively determine the consequent downstream dynamic signaling network state, which integratively governs T cell functional responses.

Our findings tangibly illustrate the crucial requirement for quantitative and dynamic integration of multiple dynamic signaling pathway activities downstream of TCR triggering for operational translation of membrane-proximal ligand affinity effects into cell functional responses such as IL-2 production. We were motivated to measure time courses by previous reports indicating the importance of ERK dynamics in cell decision processes (30) including T cell development and function (31, 32). Our modeling iterations demonstrate the information gain by inclusion of multiple time points. Furthermore, the biological significance of multiple transcription factors converging at the IL-2 promoter is reflected by the enhanced predictive power of a model consisting of multiple pathways simultaneously compared with a single molecule model. This suggests a productive avenue for exploring potentially increased efficiency of multitarget modulation of IL-2 for anti-inflammatory therapeutics and graft rejection.

Our study successfully used a multivariate linear regression model to relate signaling measurements to IL-2 production. The success of the linear regression in our study may reflect the lack of qualitative differences for downstream pathway activation by the different agonist peptides. However, previous work using PLSR has been successful in other receptor signaling pathways in which different ligands elicit binary phenotypic responses (11), suggesting that the lack of a threshold response in our system is not a key determinant for the success of linear regression modeling. It may be that the most dramatic nonlinearities will be seen in the cue-signal facet of the overall cue-signal-response paradigm, whereas the signal-response facet is more linear (15, 16, 17). Cue-signal relationships may need to be nonlinear because the upstream stimulatory cue effects typically “fan out” across multiple pathways. The signal-response relationships can be effectively approximated in terms of multilinear models because the downstream signals are integrating information across these pathways to cell phenotypic operations (17). We note that our studies involve cell population-averaged measurements, so a quantitative projection to similar network integration of receptor/ligand cues in individual cells remains uncertain; however, previous results in another system have shown that relationships of signals to responses can translate satisfactorily at least to behaviorally distinct subpopulations (12).

Although Akt is not presently known to directly activate the known IL-2 transcription factors, through our computational modeling the novel finding emerged that Akt phosphorylation levels played a leading role in accurate prediction of an IL-2 response, in context of the other key principal component signals, in a manner similar to a well-established IL-2 regulator, ERK. Prior work has noted association of costimulation-dependent IL-2 production with Akt activation (33); this may explain our observation where a strong agonist overcomes the requirement for costimulation by activating Akt, resulting in IL-2 production. Significantly, the model was able to predict a nonadditive effect between MEK/ERK and PI3K/Akt that diminished further reduction of IL-2 production by the combinatorial inhibition of both pathways. It is possible that our results in which Akt emerged as a major indicator of TCR-induced IL-2 production are due to the concerted influence of PI3K/Akt on NF-κB, NFAT, and AP-1 activation (33, 34, 35, 36, 37, 38, 39), consistent with our findings that combinatorial contributions of signaling pathways determine a cellular decision process.

We express deep appreciation to Cynthia Stokes, Arup Chakraborty, Darrell Irvine, Chris Dillon, Dinah Singer, and Al Singer for valuable input, and the anonymous reviewers for very helpful critique.

The authors have no financial conflict of interest.