Research Article

A Systems Model of Signaling Identifies a Molecular Basis Set for Cytokine-Induced Apoptosis

See allHide authors and affiliations

Science  09 Dec 2005:
Vol. 310, Issue 5754, pp. 1646-1653
DOI: 10.1126/science.1116598


Signal transduction pathways control cellular responses to stimuli, but it is unclear how molecular information is processed as a network. We constructed a systems model of 7980 intracellular signaling events that directly links measurements to 1440 response outputs associated with apoptosis. The model accurately predicted multiple time-dependent apoptotic responses induced by a combination of the death-inducing cytokine tumor necrosis factor with the prosurvival factors epidermal growth factor and insulin. By capturing the role of unsuspected autocrine circuits activated by transforming growth factor–α and interleukin-1α, the model revealed new molecular mechanisms connecting signaling to apoptosis. The model derived two groupings of intracellular signals that constitute fundamental dimensions (molecular “basis axes”) within the apoptotic signaling network. Projection along these axes captures the entire measured apoptotic network, suggesting that cell survival is determined by signaling through this canonical basis set.

Despite extensive molecular-level information on how external stimuli affect cell fate, there is minimal understanding of how such intracellular processing occurs at a systemwide level. Most extracellular “inputs” initiate complex signaling patterns that propagate through an intracellular network to change the response “outputs” that determine a cell's phenotype (1). Molecular signaling through this network is branched (fig. S1) and dynamically interconnected to the molecular history of the previous inputs, signals, and outputs (2). Thus, we sought to develop a mathematical formalism to connect signals and outputs in such a way that cellular responses could be predicted from molecular signaling patterns alone.

To reduce the biological complexity of cellular signal processing, experimental systems usually monitor changes in one defined extracellular stimulus (for example, a soluble cytokine) and one output response (such as apoptosis). However, physiological input stimuli are not processed in isolation, because signaling networks constantly receive additional inputs from changing environmental conditions, such as nutrient availability, cell density, and exposure to the extracellular matrix (3). To investigate multiinput signal processing, large-scale “systems biology” approaches have been proposed. However, most of the analyses to date have concentrated on characterizing signals down-stream of individual inputs (4, 5).

Studying multiple input stimuli requires information about the network as a whole; otherwise, intracellular changes in signal transduction molecules (called “molecular signals” hereafter) can appear paradoxical. For example, c-Jun N-terminal kinase (JNK) is a protein that has been reported to be proapoptotic (6), antiapoptotic (7), or uninvolved in apoptosis (8) in different cell systems. To investigate this JNK-apoptosis link further, we added multiple combinations of tumor necrosis factor (TNF) and epidermal growth factor (EGF) to HT-29 human colon adenocarcinoma cells and then measured the amounts of phosphorylated JNK (P-JNK) and apoptosis (Fig. 1A) (9). A plot of the TNF and EGF input stimuli and corresponding P-JNK signal and apoptosis outputs established a four-dimensional signal-response “surface” (Fig. 1B). “Slices” through this surface mimicking single TNF or EGF inputs could recapitulate any of the previously reported correlations between P-JNK and apoptosis (Fig. 1, C to E) (68). This indicated that individual molecular signals like P-JNK cannot uniquely determine a cell's commitment to apoptosis, a complex output response. Quantitative experiments that dynamically sample many critical signals would be needed (10).

Fig. 1.

JNK phosphorylation after multi-input cytokine stimulation fails to correlate with apoptotic responses. (A) JNK phosphorylation and apoptosis in HT-29 cells treated with TNF, EGF, or both. P-JNK was analyzed at 15 min by quantitative Western blotting with tubulin as a loading control, and apoptosis was measured at 24 hours by cleaved caspase-cytokeratin staining and flow cytometry (9). (B) Response surface for P-JNK (z axis) and apoptosis (color bar) induced by the input described in (A). (C to E) Lack of a correlation between JNK phosphorylation and apoptosis. P-JNK appeared to be (C) proapoptotic, (D) antiapoptotic, or (E) uninvolved in apoptosis depending on the experimental “slice” through the TNF-EGF signaling space. (F) Experimental design for TNF, EGF, and insulin stimulation. Letters correspond to the experimental treatment, used in the respective figure panels, for which 19 molecular signals (Table 1) (11) and apoptotic outputs were measured. In addition, signals and apoptotic outputs were measured in cells stimulated with 0.2-ng/ml TNF and 1-ng/ml insulin (9). The carrier was always 0.02% dimethyl sulfoxide (DMSO). (G to O) Apoptotic outputs by flow cytometry (9). In (G), (N), and (O), the insets depict representative density plots of the cleaved caspase+-cleaved cytokeratin+ population. All flow cytometry measurements were normalized to the maximum observed apoptosis for phosphatidylserine exposure (25%), membrane permeability (67%), nuclear fragmentation (21%), and cleaved caspase-cytokeratin (77%). Data are presented as the mean ± SEM of triplicate biological samples in panels (C) to (E) and (G) to (O).

Dense experimental sampling of 7980 molecular signals and 1440 response outputs for cytokine-induced apoptosis. To investigate how signaling networks coordinate cellular output responses, we used an established multi-input system where HT-29 cells are stimulated with three biologically relevant cytokines: TNF, EGF, and insulin (11, 12). The intracellular protein network downstream of these cytokine inputs is understood in reasonable detail (fig. S1). Within well-recognized network branches, 19 intracellular measurements of key receptor, kinase, caspase, and adaptor proteins were collected (Table 1 and fig. S1) (11). Using nine distinct pairwise combinations of TNF, EGF, and insulin (Fig. 1F), we sampled each molecular signal in triplicate at 13 time points between 0 and 24 hours to compile 7980 distinct molecular signals from the shared intracellular network (11). TNF and EGF (or insulin) are opposing cytokine inputs that promote or inhibit apoptosis, respectively (13, 14), and the intracellular measurements thus identified which molecular signals were activated en route to the cellular decision to die or to survive. However, it remained unclear how EGF- and insulin-induced molecular signals antagonized TNF-induced apoptosis.

Table 1.

Signaling network measurements. Biological functions were assigned based on the most-recognized property of the protein. Molecular signals indicate the biochemical property of the protein that was measured. Assay indicates whether the protein was measured by high-throughput kinase activity assay (10), antibody microarray (33), or quantitative Western blotting (12).

Protein name Biological function Molecular signal Assay
IKK Ser kinase Kinase activity Kinase assay
JNK1 Ser-Thr kinase Kinase activity Kinase assay
MK2 Ser-Thr kinase Kinase activity Kinase assay
EGFR Receptor Tyr kinase Phosphorylation (Tyr1068) Ab microarray
Total amount Ab microarray
Phospho/total ratio Ab microarray
MEK Dual-specificity kinase Phosphorylation (Ser217/Ser221) Western blot
ERK Ser-Thr kinase Kinase activity Kinase assay
IRS1 Adaptor-scaffold Phosphorylation (Ser636) Western blot
Phosphorylation (Tyr896) Western blot
Akt Ser-Thr kinase Phosphorylation (Ser473) Ab microarray
Total amount Ab microarray
Kinase activity Kinase assay
Phosphorylation (Ser473) Western blot
Phospho/total ratio Ab microarray
FKHR Transcription factor Phosphorylation (Ser256) Western blot
Caspase-8 Cys protease Zymogen amount Western blot
Cleaved amount Western blot
Caspase-3 Cys protease Zymogen amount Western blot

To test whether apoptosis could be connected to the measured signaling network (11), we measured the cell-death phenotype for each combinatorial cytokine stimulus (Fig. 1F). The apoptotic response itself involves various cellular changes (outputs) that can be regulated independently (15). Individual parameters used to characterize cell death (such as loss of membrane asymmetry) only partially reflect the overall cellular response. Therefore, we selected four distinct apoptotic outputs (phosphatidylserine exposure, membrane permeability, nuclear fragmentation, and caspase substrate cleavage) and measured each output response by flow cytometry 12, 24, and 48 hours after stimulation (Fig. 1, G to O, and database S1) (9). Together, these output measurements constituted an apoptotic “signature” that characterized early (phosphatidylserine exposure), middle (caspase substrate cleavage and membrane permeability), and late (nuclear fragmentation) responses of apoptosis (the biological outputs in our system).

Apoptotic signatures were measured for the nine cytokine input combinations (Fig. 1F), revealing temporal and cytokine dose-dependent features that would have been missed by examining single apoptotic outputs alone. Membrane permeability (red) and caspase substrate cleavage (purple) measured dead cells cumulatively (16). Therefore, these outputs increased monotonically with time and TNF dose (Fig. 1, G, J, and M). In contrast, phosphatidylserine exposure (loss of membrane asymmetry before membrane permeability) and nuclear fragmentation (digestion of DNA before complete cellular fragmentation) were transient cell states (16) that could increase and decrease with time and input dose (Fig. 1, G to O).

High concentrations of EGF and insulin antagonized TNF-induced cell death in a similar manner, particularly reducing apoptotic responses at early times (Fig. 1, M to O, and fig. S2). Low concentrations of EGF and insulin elicited different output responses. EGF reduced TNF-induced membrane permeability at 48 hours but increased transient phosphatidylserine exposure at 12 hours (Fig. 1, J and L; P < 0.05, Student's t test with Bonferroni correction). In contrast, insulin reduced TNF-induced membrane permeability at 12 and 48 hours, and phosphatidylserine exposure was decreased at 12 hours but higher at 48 hours compared with TNF alone (Fig. 1, J and K; P < 0.05). These measured apoptotic signatures provide evidence that different apoptotic outputs can be controlled separately, depending on the input stimulus. For example, membrane permeability is most commonly associated with secondary production of reactive oxygen species (17), whereas phosphatidylserine exposure is thought to involve caspase- and Ca2+-dependent processes (18). The JNK family protein kinases (Fig. 1A) have been implicated in reactive oxygen signaling (19, 20) but are thought to be largely independent of Ca2+ signaling (21). Thus, relating the intracellular network dynamics to the complete apoptotic signature requires more information than can be provided by individual molecular signals, like JNK1.

A partial least-squares model predicts the 12 cytokine-induced apoptotic outputs from a 660-element signaling vector. We sought to determine whether multiple molecular signals in combination could quantitatively capture the entire apoptotic signature as a response and predict apoptosis more globally. We designed a mathematical formalism that could identify the information content within each molecular signal that most closely mapped onto the output responses. The resulting mapping of lumped signals to corresponding responses might then allow us to identify the most relevant “information variables” for apoptosis and use these variables to predict apoptotic responses to stimuli outside the training set.

It was unclear what aspect of a dynamically sampled molecular signal should fill the role of an information variable for apoptosis. For example, in addition to measurements of kinase activity at various time points, it was not known if other data such as the maximum activity, the rate of rise of activity, or the time when peak activity occurred also contained useful information. Therefore, we defined a panel of time-dependent signaling “metrics” that could be derived empirically from any dynamic signal (Fig. 2A). Each of the metrics [for example, activity of the mitogen-associated protein kinase (MAPK) extracellular signal–regulated kinase (ERK) at 5 min, ERK activity at 15 min, peak ERK activity, and the abundance of the phosphorylated (active) MAPK kinase MEK (P-MEK) at 5 min] can then be represented by an “axis” along which particular multi-input stimuli project (Fig. 2C).

Fig. 2.

A data-driven principal components-based model correctly predicts apoptosis and caspase activation from molecular signals activated by multiinput TNF, EGF, and insulin stimulation. (A) Extraction of signaling and response metrics. (Top) Time-course measurements (11) of the signaling network were used to group individual time points (t), instantaneous derivatives (d/dt), area under the curve (AUC), and the maximum (Max), mean, and steady-state (S-S) values and to form a signaling metric set, in this example for ERK. The activation slope and decay rates for signaling peaks were also incorporated (9). (Bottom) Individual apoptotic output time points were concatenated for each treatment into an apoptosis vector. (B) Construction of the PLS apoptosis model. Individual signaling metric sets were concatenated into a single signaling vector for each treatment condition, and these vectors then regressed against the corresponding apoptotic response vectors (22). wA is the coefficient matrix of the apoptosis regression model. (C) PLS-based supervised decomposition of signaling vector space into principal components. (Top) A simplified example of dimensionality reduction by principal components analysis. In this example, collinear axes such as peak ERK activity, ERK activity at 5 min, and P-MEK levels at 5 min could be reduced into a single principal component axis (p1). (Bottom) Principal component axes were captured in a supervised manner to maximize prediction of apoptotic responses. In the full PLS model, the four apoptotic outputs (red, green, blue, and purple quandrants of each data point) at three time points (12, 24, and 48 hours slices within quadrants) were oriented along three principal component axes (p1, p2, and p3) determined to contain the essential signaling information for predicting apoptosis. Each point corresponds to one multi-input stimulus, such as TNF+EGF, and the shading of the point corresponds qualitatively to the extent of measured apoptosis. The numbered axes (left) indicate separate metrics extracted from different molecular signals, shown in (A). (D) RMSE of calibration and prediction of apoptotic outputs by the PLS model as a function of increasing number of principal components. An optimum model with three components was selected. (E to G) Correlation plots between measured phosphatidylserine exposure (y axis) and cross-validated predictions of phosphatidylserine exposure (x axis) by the PLS model including (red) and not including (yellow) caspase signals. The merged overlay is shown in green. Marker size corresponds to the response time point at 12 (Embedded Image), 24 (Embedded Image), and 48 (⚫) hours. Data are presented as the mean ± SEM, and model uncertainties were estimated by jack-knifing (22). (H to J) Correlation plots between measured and cross-validated predictions of membrane permeability, nuclear fragmentation, and caspase-substrate cleavage. Data are presented as in (E) to (G). (K) Comparisons between measured caspase levels and cross-validated predictions from a specific PLS model of caspase activation. Signaling metric vectors from noncaspase molecular signals were regressed against the 0- to 24-hour time-point measurements of procaspase-8, cleaved caspase-8, and procaspase-3. wA is the coefficient matrix of the caspase regression model. The nine rows in the measured and predicted color coding correspond to the nine treatment combinations shown in Fig. 1F. Each region of the heat map was normalized separately for comparison.

In total, we extracted 30 to 40 metrics (table S1) from each time course of a molecular signal to form a composite “metric set” (Fig. 2A). Metric sets were defined for all 19 molecular signals, forming a group of individual axes that together define a 660-dimensional signaling space. The projection of a multi-input stimulus along these axes creates a “signaling vector” of 660 column elements corresponding to the 30 to 40 metrics from each of the 19 measured molecular signals (Fig. 2B). The apoptotic outputs were similarly concatenated to form an apoptosis vector, where the 12 column elements are the four apoptotic outputs at three time points each (Fig. 2A). Each input stimulus has its own particular signaling and apoptosis vectors, so these were calculated separately for each of the nine treatment combinations (Fig. 1F).

We sought to identify the best mapping of the 660-dimensional signaling metric space onto the 12-dimensional output response space. A simple linear mapping would require the 7920-element coefficient matrix (12 rows by 660 columns) that best transformed the signaling space into the response space (Fig. 2B), a calculation that would be impossible to solve uniquely given only nine multi-input conditions. To simplify the signal-to-response mapping and retain biological meaning, we first assumed that some of the signaling metric dimensions were either redundant or irrelevant. For instance, the peak activity of ERK contains the same information as the ERK activity at 5 min, when the ERK peak often occurred (Fig. 2C). Likewise, other axes (e.g., metrics at 0 min, before the stimulus) merely scrambled the projections because they pointed in directions unrelated to apoptosis (Fig. 2C). These redundant and uninformative axes could be deemphasized or eliminated without any loss of information about the apoptotic outputs. Second, we assumed that the remaining informative axes could be compressed by linear combination into a small number of dimensions that retained the critical apoptotic information and thus constituted a biologically relevant basis set for the signal-to-response mapping. For instance, proteins like MEK and ERK are part of the same signaling pathway and were thus activated similarly. Although not identical, the information in many MEK and ERK metrics pointed in collinear directions that could be combined into a MEK-ERK “super axis” that retained the projection with fewer dimensions (Fig. 2C).

To reduce unnecessary axes and condense important axes mathematically, we used partial least-squares (PLS) regression, which simplifies dimensions on the basis of their covariance with a specified dependent variable (22, 23). The original 660-dimensional signaling space was reduced to a series of super axes that together best orient the measured apoptotic-output elements within the apoptosis vector (Fig. 2, A and C). PLS modeling, like singular value decomposition (24), calculates super axes as an orthogonal set of “principal components,” which contain linear combinations of the original 660 metric dimensions weighted by their contribution to the apoptotic outputs (Fig. 2C). Principal components are calculated iteratively so that successive PLS dimensions are regressed against the apoptotic-output information not captured by the preceding component. After several iterations, including more dimensions becomes undesirable, because the residual information is so small that new principal components capture spurious fluctuations in the outputs, such as measurement error and noise, rather than meaningful data (22).

To optimize the number of model dimensions, we examined the root-mean-squared error (RMSE) between the measured apoptosis vector and the values from models with increasing numbers of principal components. All of the input treatments were included in the initial model training to assess the RMSE of data fitting (calibration). Each treatment was then individually withheld from the training set to construct a cross-validation model, in which an RMSE of prediction could be assessed by predicting the withheld sample. The calibrated RMSE decreased monotonically, but the predicted RMSE was minimized with just three principal components (Fig. 2D). Using the resulting three-component apoptosis model, we examined the correlation between the measured apoptotic outputs and the cross-validated predictions for each treatment (22). If the correct signaling information was retained by the principal components in the model, then the apoptotic outputs for TNF, EGF, and insulin stimuli not included during the model training should also be predicted. Indeed, we found a high correlation for all 12 apoptotic outputs, with predictions that were accurate to within 94% of the measured values overall (Fig. 2, E to J). Although the 12 individual apoptotic outputs differed quantitatively from one another (Fig. 1, G to O), their stimulus-dependent changes were almost entirely captured by a single PLS model of the intracellular network (Fig. 2, E to J). The 660 dimensions of the original signaling space (Fig. 2C) had thus been condensed to three dimensions, which were enriched in the molecular signals most useful for predicting the apoptotic outputs. Furthermore, specifying the dynamic state of the intracellular signaling network was absolutely essential—an equivalent PLS model given only cytokine-input concentrations (Fig. 1F) rather than intracellular signaling metrics could not correctly predict the apoptotic outputs (45% accuracy) (25).

The original measurements (11) of the signaling network contained several direct effectors of apoptosis, such as caspases (Table 1) (26). The model would obviously be less valuable if these late-effector signals were providing all of the predictive power to the model. To test whether caspase metrics were required for accurate predictions, we removed all of the caspase signals from the model and rederived the principal component axes. The resulting predictions were essentially identical (Fig. 2, E to J). Thus, non-caspase signals in the network contained more than enough information to predict the apoptotic outputs quantitatively. Furthermore, molecular signals activated before the initial onset of apoptosis were themselves sufficient for predicting the apoptotic signature. We found that a separate PLS apoptosis model derived exclusively from signaling measurements made at 0 to 4 hours after cytokine addition was accurate to within 81% (fig. S3).

To examine the relationship between other molecular signals and late effector caspases directly, we next removed the apoptotic outputs altogether and defined the procaspase-8, cleaved caspase-8, and procaspase-3 time-point measurements as a new set of cellular outputs (Fig. 2K). Using the remaining network measurements, the PLS model predicted the caspase response dynamics within 81% accuracy (Fig. 2K). Together, these results suggested that both the caspase effector signals and the final cellular output responses were encoded by the upstream signaling network.

Model-driven discovery of apoptosis regulation by means of autocrine circuits. In HT-29 cells, two regulated autocrine stimuli—transforming growth factor–α (TGF-α) and interleukin-1α (IL-1α)—cooperate with TNF to activate molecular signals in the network (Fig. 3A) (11). Whether these autocrine circuits contribute substantially to TNF-induced apoptosis is not known. We therefore used the PLS model to predict what the apoptotic signature would be when autocrine TGF-α and IL-1α circuits were disrupted with either an antibody to the EGF receptor (C225) or an IL-1 receptor antagonist (IL-1ra), and the cells were then treated with TNF. All 19 molecular signals were measured from 0 to 24 hours (11) and provided as input data to the PLS model. We then measured the actual TNF+C225- and TNF+IL-1ra–induced apoptotic signatures experimentally (9) and compared these with the model predictions. This experiment was a particularly stringent test of the model, because neither TGF-α nor IL-1α had been part of the original training set (Fig. 1F).

Fig. 3.

The principal components-based model captures hidden autocrine feedback in the signaling network and identifies late MK2 activity as a TGF-α–induced prosurvival signal. (A) Diagram of TNF-induced autocrine circuits in HT-29 cells (11). Orange arrows indicate fast pathways (before 4 hours) and purple arrows indicate slow pathways (after 4 hours). C225 and IL-1ra were used as pharmacological inhibitors of autocrine TGF-α and IL-1α, respectively. (B) The model correctly predicts TNF-induced apoptosis when the autocrine feedback circuits are disrupted by C225 and IL-1ra. Apoptotic outputs were measured as in Fig. 1, G to O, for cells treated with 5-ng/ml TNF and 10-μg/ml C225 (circle) and 100-ng/ml TNF plus 10-μg/ml IL-1ra (square). Marker color corresponds to the apoptotic index (as in Fig. 1, G to O) and size corresponds to the time point (as in Fig. 2, E to J). Horizontal error bars indicate model uncertainty by jack-knifing (22). (C) The model correctly predicts caspase cleavage in response to TNF stimulation when the autocrine feedback circuits are disrupted. Comparisons between measured caspase levels and model predictions for the autocrine circuit perturbations. Each panel was normalized separately, and the unperturbed caspase dynamics were included for comparison. (D) The TNF→IL-1α autocrine circuit is proapoptotic, whereas the TNF→TGF-α circuit sends both pro- and antiapoptotic signals. Comparison of TNF-induced apoptotic outputs with and without autocrine IL-1α and autocrine TGF-α. Data are shown as percent change in apoptosis in C225- and IL-1ra–treated cells compared with unperturbed TNF-stimulated cells (10). (E) TNF-induced late MK2 activity is TGF-α dependent. Cells were stimulated with 5-ng/ml TNF in the presence or absence of 10-μg/ml C225 to block autocrine TGF-α, and MK2 activity was measured by in vitro kinase assay (10). (F) Late MK2 activity is restricted to surviving cells. (Left) Cells were treated with TNF for 20 hours, separated into floating and adherent subpopulations, and MK2 activity was assayed as in (E). Data were normalized to basal MK2 activity as described (12). (Right) Floating and adherent cells correspond to apoptotic and viable cells, respectively, as indicated by caspase-8 cleavage and 4′,6′-diamidino-2-phenylindole (DAPI) staining (25). (G) Late MK2 activity is a prosurvival signal. Cells were treated for 24 hours with 100 ng/ml TNF in the presence or absence of the inhibitor SB202190 (SB) during the indicated times and assayed for apoptosis by the presence of caspase substrate cleavage. For (B) and (D) to (G), data are presented as the mean ± SEM of triplicate biological experiments (11).

We observed a 90% correlation between the measured apoptotic outputs and the model predictions when the two autocrine loops were disrupted individually (Fig. 3B), indicating that the model could predict the contributions of cytokines other than TNF, EGF, and insulin. Furthermore, upstream molecular signals alone were highly predictive of downstream effector caspase activation after TNF+C225 and TNF+IL-1ra stimuli (prediction within 91% of measured values; Fig. 3C). Therefore, the contributions of autocrine TGF-α and IL-1α to TNF-induced apoptosis had been implicitly and correctly incorporated throughout models of apoptotic outputs as well as caspase activation.

Disruption of the IL-1α autocrine loop decreased most apoptotic outputs in cells responding to TNF (Fig. 3D), suggesting that autocrine IL-1α was an extracellular positive-feedback circuit for TNF-induced apoptosis. In contrast, disruption of TGF-α signaling by C225 led to large changes in TNF-induced activation of the network but did not lead to any clear overall changes in apoptosis (Fig. 3D). Some outputs, such as phosphatidylserine exposure, were increased in cells treated with C225 (P < 10–5, two-way analysis of variance). Others, such as caspase substrate cleavage, were decreased (P < 0.05), whereas 6 of the 12 individual apoptotic outputs did not change in a statistically significant manner (9, 25).

TGF-α is a member of the EGF family (27), and treatment of cells with EGF was shown to reduce apoptosis (Fig. 1, M and O). Why then did blocking autocrine TGF-α with C225 not increase apoptosis overall? First, in HT-29 cells, autocrine TGF-α activates a prodeath signal by inducing late (12- to 24-hour) release of autocrine IL-1α (Fig. 3A and fig. S4) (11). Second, we reasoned that the late, prodeath IL-1α signal must be offset by some unknown late, prosurvival signal. Concurrent prodeath and prosurvival signals at late times could then neutralize one another and explain the net lack of an effect of autocrine TGF-α on apoptosis. To help identify this late, prosurvival signal, we examined the coefficients within the PLS model that quantified the contribution of each signal to apoptosis (Fig. 2B), because the model correctly predicted the mixed apoptotic response resulting from autocrine TGF-α perturbation (Fig. 3B). Early MAPK-activated protein kinase 2 (MK2) kinase activity in the model correlates with apoptosis (25); by contrast, the model revealed that MK2 signaling at times after 12 hours was strongly anticorrelated with apoptosis, implicating MK2 as a prosurvival kinase at later times. In agreement with this, when autocrine TGF-α was blocked by C225, MK2 activity at late times was inhibited (Fig. 3E) (11, 25). Thus, both the model and experiment suggested that late MK2 activity was a prosurvival signal activated by autocrine TGF-α.

If effectively prosurvival, then MK2 activity at late times should be present primarily in viable cells. MK2 signaling was therefore analyzed separately in live and apoptotic cells by measuring the floating and adherent subpopulations of TNF-treated cells. All MK2 activity was detectable only in the adherent (viable) cells (Fig. 3F). If we permitted the early phase of MK2 activity to occur in response to TNF but rapidly inhibited the late phase with SB202190, a small-molecule inhibitor of the MK2-activating kinase p38 (28, 29), we observed a significant increase in cell death at 24 hours (Fig. 3G; P < 10–5, Students t test). Taken together, these experiments demonstrate that the signal-to-response linkages implicated by the PLS model can reveal new biological mechanisms that would not be easily recognized without a mathematical formalism.

Redundant encoding of signaling components and critical roles of MAPKs in PLS predictions of apoptosis. In the full model, 660 metrics derived from 7980 signaling measurements were used to predict apoptosis (Fig. 2, E to J, and Fig. 3B). In retrospect, it was not clear whether apoptotic outputs could be predicted equally well by a reduced number of metrics derived from smaller, more tractable experiments. We found that a model containing only the 20 most informative metrics, as determined by the relative magnitude of their coefficients in the model (table S2) (9), was nearly as predictive of apoptosis as one that used all 660 metrics (Fig. 4A). A noteworthy feature of the top 20 metrics was that they were not the obvious metrics we would have chosen on the basis of our basic understanding of the regulation of apoptosis. Activation slopes of Akt and insulin receptor substrate 1 (IRS1) phosphorylation and integrated peaks of JNK1 and inhibitor of nuclear factor κB kinase (IKK) signaling were included, but caspase metrics were missing entirely (table S2). Use of fewer than these top 20 metrics, in isolation, gave substantially less effective predictions (25), supporting our original hypothesis that individual signaling measurements would not broadly predict cell responses (Fig. 1, C to E).

Fig. 4.

Two orthogonal stress and survival axes link complex stimuli to apoptosis. (A and B) The principal component-based model is redundantly encoded. (A) Decreased accuracy of prediction of apoptotic responses to TNF, EGF, and insulin with a decreasing number of signaling metrics. (B) Decrease in accuracy of prediction of apoptotic responses to autocrine signals with a decreasing number of molecular signals. Both metrics and molecular signals were eliminated sequentially from best to worst on the basis of their variable importance in the principal-component projection (9). (Inset) Predictive efficiency as a function of number of molecular signals included. To calculate efficiency, the TNF-EGF-insulin apoptotic predictions in (B) were divided by the number of molecular signals in the submodel. For (A) and (B), data are presented as the central prediction ± 90% Fisher Z-transformed confidence intervals. (C) The two principal components can correspond to distinct time-resolved mechanisms of molecular-signal activation. Contribution of IKK activity time-point metrics to axes 1 and 2. Early and late IKK time points are shaded in light and dark gray, respectively. A reciprocal time-dependent change in projection along axes 1 and 2 was observed for MK2 (25). (D to F) The intracellular signaling network is defined by two orthogonal stress and survival axes that link complex stimuli to apoptosis. (D) Cytokine inputs project along orthogonal stress and survival axes (axes 1 and 2) of the PLS model. (E) Cytokine combinations project differently from linear superposition of isolated cytokine stimuli. EGF in combination with TNF directly antagonizes stress pathway signaling, whereas insulin in combination with TNF both antagonizes stress signaling and induces separate prosurvival signaling pathways. The linear superposition of single-input TNF, EGF, and insulin projections from Fig. 4E is shown in gray. (F) Autocrine feedback circuits rotate the TNF-induced signaling network to reinforce stress responses and suppress any signaling down the survival axis. For (D) to (F), samples were projected using the model scores as described (22).

High predictive ability of the top 20 metrics (Fig. 4A) suggested redundancy in the signaling information contained within the original 660-metric model. To investigate this, we sequentially removed various portions of the most informative metrics and then recalculated the prediction of apoptosis by cross-validation. Up to the top 350 metrics could be eliminated from the full list of 660 metrics before the apoptosis model lost all predictive ability (Fig. 4A). This implied that the full set of 660 biological metrics was redundantly encoded with the stimulus-specific information required to mediate all of the apoptotic outputs (30, 31). We also analyzed the contributions of individual proteins (that is, the 19 molecular signals) based on the average information contained in their derived metrics. The top three molecular signals—JNK1, MK2, and ERK—belonged to MAPK cascades (25). A model given only metrics from these three signals performed nearly as well as the full PLS model in predicting the effects of autocrine perturbations (Fig. 4B). A model derived from the seven least informative molecular signals was also predictive (Fig. 4B). Together, these findings suggest that measurements of three to seven relevant intracellular molecular signals are sufficient to predict network-dependent output responses. Prediction efficiency (defined as the relative predictive ability divided by the number of molecular signals included in the model) was maximal with four to five molecular signals (Fig. 4B, inset).

The PLS model identifies stress-apoptosis and survival signaling axes. We examined the relationship between all of the signaling network measurements and the input treatments as signals propagated through the model's first two principal components (22), which predicted the apoptotic outputs within 92%. These components are the two main lumped contributions of the signaling metrics, which form a pair of orthogonal axes defining the optimal two-dimensional slice through the signaling data set (22). We found that certain signals and treatments were clearly overrepresented in these dimensions. The first principal component, axis 1, was oriented toward stress and apoptotic pathways and included early JNK1 activity, early MK2 activity, and late cleaved caspase-8 metrics (table S3). In contrast, the second principal component, axis 2, appeared to constitute a global survival signal that included phosphorylated Akt (P-Akt), phosphorylated IRS1 (P-IRS1), phosphorylated Forkhead transcription factor (P-FKHR), and procaspase-3 metrics (table S4). The specific molecular signals emphasized in these components were entirely consistent with the known molecular mechanisms that link these signaling proteins to apoptosis (fig. S1).

Using the principal component axes, we could reanalyze the signaling contributions to apoptosis from cells exposed to single or combined cytokines. Certain signals, such as IKK, contributed differently to the stress-apoptosis and survival axes depending upon the time point when the molecular signal was measured. Early IKK activity induced directly by the TNF receptor (TNFR) was weighted predominantly along prosurvival axis 2 (Fig. 4C). However, IKK activity after 12 hours, which occurs indirectly in response to an autocrine signal from IL-1α (Fig. 3A) (11), contributed more to the stress-apoptosis axis (Fig. 4C). Thus, the PLS model had revealed that the same molecule, such as IKK, can convey either pro- or antiapoptotic messages, depending on the timing and mechanism of activation.

For individual cytokines, we found that TNF treatment alone projected strongly along prodeath axis 1, whereas isolated EGF and insulin treatments mapped exclusively on prosurvival axis 2 (Fig. 4D). This reinforced our original intuition that TNF and EGF-insulin stimuli act upon orthogonal and antagonistic signaling axes for apoptosis. In contrast, the multi-input projections were markedly different from what would be predicted by a summation of the single-input treatments (Fig. 4E, gray). TNF, EGF, and insulin each lost a fraction of their original projections along the two axes, indicating that the input stimuli were antagonized when added in combination (Fig. 4E). The TNF+EGF and TNF+insulin treatments were distinctly separated from one another: EGF appeared to antagonize TNF-induced apoptosis by specifically reducing the projection along the stress-apoptosis axis 1 without any change along axis 2 (Fig. 4E); in contrast, insulin actively promoted prosurvival signaling along axis 2 while also inhibiting stress-apoptosis signaling along axis 1. Therefore, analyzing the multi-input stimuli through these model-derived biological “basis axes” (Fig. 4D) helped to reveal the different network strategies used by EGF and insulin to antagonize TNF-induced apoptosis.

Finally, to determine the contributions of TNF-induced autocrine circuits in the model, we mapped the TNF+C225 and TNF+IL-1ra treatments (Fig. 4F). The projection of TNF along the stress-apoptosis axis (Fig. 4F) was enforced by the autocrine circuits, which increased the contribution along axis 1 and decreased the contribution along axis 2. This is consistent with the notion that regulated autocrine circuits provide microenvironment-dependent feedback to cells during phenotypic decision processes, such as death-survival (11, 32). Furthermore, it illustrated directly that effects of complex environmental stimuli were entirely contained within the two canonical basis axes distilled from the original 660-dimensional signaling metric space by the PLS model (Fig. 4D).

In summary, by using a systems approach that combines quantitative experiments with data-driven modeling, we identified two canonical axes—a stress-apoptosis axis and a survival axis—that together constitute a molecular basis set for the signaling network that controls apoptosis. These axes capture the dynamic intracellular signal processing of diverse stimuli, including autocrine-feedback circuits. Our work illustrates how a complex signaling network can be reduced empirically to a much simpler computational model that is directly tied to biological mechanism.

Supporting Online Material

Materials and Methods

Figs. S1 to S4

Tables S1 to S4

Database S1


References and Notes

View Abstract

Stay Connected to Science

Navigate This Article