Research ArticleCancer

Signaling pathway models as biomarkers: Patient-specific simulations of JNK activity predict the survival of neuroblastoma patients

See allHide authors and affiliations

Science Signaling  22 Dec 2015:
Vol. 8, Issue 408, pp. ra130
DOI: 10.1126/scisignal.aab0990

Pathway activity as a biomarker

Understanding signaling networks may enable prediction of disease mechanisms or responses to therapeutic strategies. Fey et al. showed that constructing a pathway that reproduces the all-or-nothing, switch-like activation of the stress-activated kinase JNK could be used to stratify neuroblastoma patients. Switch-like activation of JNK leads to cell death. By integrating patient-specific information about the abundance of the components of the JNK pathway in neuroblastoma samples into the model, the authors simulated the activity of the pathway and accurately predicted survival on the basis of the dynamic properties of the pathway. Thus, pathway activity served as a biomarker. The results also showed that alterations in the network that prevent the switch-like activation of JNK were associated with poor survival of neuroblastoma patients, thus providing potential molecular mechanisms for the inherent resistance of some of these tumors to treatment.


Signaling pathways control cell fate decisions that ultimately determine the behavior of cancer cells. Therefore, the dynamics of pathway activity may contain prognostically relevant information different from that contained in the static nature of other types of biomarkers. To investigate this hypothesis, we characterized the network that regulated stress signaling by the c-Jun N-terminal kinase (JNK) pathway in neuroblastoma cells. We generated an experimentally calibrated and validated computational model of this network and used the model to extract prognostic information from neuroblastoma patient–specific simulations of JNK activation. Switch-like JNK activation mediates cell death by apoptosis. An inability to initiate switch-like JNK activation in the simulations was significantly associated with poor overall survival for patients with neuroblastoma with or without MYCN amplification, indicating that patient-specific simulations of JNK activation could stratify patients. Furthermore, our analysis demonstrated that extracting information about a signaling pathway to develop a prognostically useful model requires understanding of not only components and disease-associated changes in the abundance or activity of the components but also how those changes affect pathway dynamics.


Traditional cancer biomarkers are based on an association between clinical parameters and the expression or mutation of specific genes, often combined into panels to strengthen their sensitivity and specificity. Although these biomarkers have extensive clinical utility, they frequently lack any mechanistic anchoring to the fundamental cellular processes responsible for tumorigenesis or therapeutic response. Signaling pathways play a key role in these processes; therefore, analysis of the activity of signaling pathways should aid in the development of biomarkers linked to cellular function. However, it is technically challenging to observe the activation of these pathways during tumor formation or therapeutic treatment. We therefore hypothesized that simulating the potential of cancer cells to activate a central signaling pathway under these conditions may provide insight into the prognostic information contained in these pathways.

Focusing on the c-Jun N-terminal kinase (JNK) pathway, because this pathway mediates apoptotic cell death induced by both physiological stress and therapeutic agents, we sought to develop a model of the JNK pathway that could correlate patient-specific JNK response dynamics with survival data. Although the identities of the kinases that activate JNK are well known (1), this signaling pathway displays dynamic behavior suggestive of a complex network structure that belies a simple linear kinase cascade. Systems-level approaches have previously revealed a complex network of genes that influence basal JNK activity in migratory Drosophila cells (2), but the biochemical mechanisms that govern apoptotic JNK activation are still lacking. Here, we have delineated a network structure that regulates stress-induced JNK activation in the childhood tumor neuroblastoma, which undergoes JNK-dependent apoptosis in response to various stimuli (310).

The JNK pathway has a classical three-tiered mitogen-activated protein kinase (MAPK) structure and, like all MAPK pathways, mediates diverse cellular responses. The kinases in the MAPK cascade that activate JNK are MKK4 and MKK7, which in turn are activated by various MAPK kinase kinases (MAPKKKs) (1, 11). Growth factor–induced transient activation of JNK promotes cell survival and proliferation. This transient activation of JNK exhibits a graded response to stimulus (growth factor) concentration (12, 13). In contrast, stress-induced prolonged JNK activity that results in apoptosis exhibits an ultrasensitive switch-like response (1215). Ultrasensitivity can result from positive feedback, and we previously proposed a positive feedback loop between JNK and its upstream kinases as a mechanism to achieve this ultrasensitive JNK activation (16).

Neuroblastoma originates from neural crest cells, typically arising in the paravertebral sympathetic ganglia and adrenal medulla of young children. Neuroblastoma patients are routinely classified into low-, intermediate-, and high-risk groups according to disease stage, patient age, and MYCN amplification status (1719). Although considerable progress has been made in improving survival rates for intermediate-risk patients, the overall survival rate for high-risk neuroblastoma patients remains less than 40 to 50% (18, 19). Amplification of the transcription factor MYCN occurs frequently in high-risk neuroblastomas and drives aggressive tumor behavior, resistance to chemotherapy (chemoresistance), and poor patient outcome (1921). However, many patients lacking MYCN amplification also have a poor prognosis, and there is little understanding of the biochemical pathways that drive tumorigenesis and chemoresistance in these patients (22). The JNK network structure in neuroblastoma is undefined.

Here, we developed a dynamic mathematical model of the JNK network on the basis of experimental evaluation of stress-induced JNK signaling in SH-SY5Y neuroblastoma cells and used the model to generate patient-specific simulations of stress-induced JNK activation. The results revealed a central role for JNK signaling dynamics in neuroblastoma and represent a model-based biomarker that robustly stratified neuroblastoma patients across different individual molecular backgrounds.


Ultrasensitive JNK activation and apoptosis in neuroblastoma cells exposed to stress

Consistent with a tumor-suppressive role for JNK signaling in response to stress within the tumor microenvironment, we observed that JNK inhibition increased tumor growth by about threefold in a zebrafish model of MYCN-driven, spontaneous neuroblastoma (23) (Fig. 1A). This proapoptotic role of JNK also occurred in the SH-SY5Y neuroblastoma cell line, in which a JNK inhibitor dampened the apoptotic response to various stress conditions (Fig. 1B) and significantly reduced stress-dependent p53 accumulation and caspase-3 cleavage (fig. S1A) triggered by physiological or drug-induced stress. We observed this reduction in apoptotic response in SH-SY5Y cells using agents that inhibit protein translation (anisomycin) to produce ribotoxic stress, inhibit mitochondrial function (rotenone) to produce oxidative stress, cause DNA damage (doxorubicin) to produce genotoxic stress, or disrupt microtubules (vincristine) to cause cell cycle arrest (Fig. 1B and fig. S1A). Two of these are drugs used for neuroblastoma treatment: Doxorubicin is used to treat intermediate-risk neuroblastoma patients, and vincristine is used in combination therapy for high-risk neuroblastoma patients (19).

Fig. 1 Ultrasensitive, apoptotic JNK signaling in neuroblastoma.

(A) Analysis of neuroblastoma growth in zebrafish, expressing MYCN-EGFP from the dopamine β-hydroxylase promoter, in the presence of JNK inhibitor II (5 μM, 14 days) or dimethyl sulfoxide (DMSO) control (mean ± SEM, n = 4). **P < 0.01; *P < 0.05, two-sample Student’s t test). (B) Apoptosis of SH-SY5Y cells (sub-G1 population) after treatment with anisomycin (300 nM), rotenone (1000 nM), vincristine (100 nM), or doxorubicin (100 nM) for 24 hours, in the presence or absence of JNK inhibitor VIII (10 μM) (mean ± SD, n = 3). (C) Activation loop phosphorylation of JNK, p38, and ERK after treatment of SH-SY5Y cells with increasing concentrations of anisomycin (3 to 1000 nM for 30 min) (mean ± SD, n = 3). H, Hill exponent; CI, confidence interval. (D) Activation loop phosphorylation of JNK, MKK4, and MKK7 after treatment of SH-SY5Y cells as in (C) except in the presence or absence of JNK inhibitor II (20 μM) (mean ± SD, n = 3). **P < 0.01; *P < 0.05, two-sample Student’s t test. (E) Activation loop phosphorylation of JNK, MKK4, and MKK7 after treatment of SH-SY5Y cells with control or ZAK siRNA for 48 hours, followed by treatment with increasing doses of the indicated stimuli. Data are representative of three experiments. In panels (C) and (D) quantified data, phosphorylation is calculated relative to the maximal observed value for each phosphosite.

To determine whether the apoptotic, ultrasensitive mode of JNK activation occurs in neuroblastoma cells, we exposed SH-SY5Y cells to increasing concentrations of anisomycin for 30 min (Fig. 1C), which we found, was the time point representing the maximal abundance of the phosphorylated and active form of JNK (fig. S1, B and C). The antibody used recognizes all three JNK isoforms JNK1, JNK2, and JNK3 (also called MAPK8, MAPK9, and MAPK10) phosphorylated in the activation loop (pJNKT184/Y185). Anisomycin induced ultrasensitive JNK activation, with pJNKT183/Y185 increasing from basal to saturated amounts within one order of magnitude change in anisomycin concentration. We explored the response of two other MAPKs to anisomycin, the stress-activated kinase p38 and extracellular signal–regulated kinase 1 (ERK1) and ERK2 (collectively referred to as ERK). In contrast to the change in the abundance of pJNKT183/Y185, the abundance of activation loop–phosphorylated p38 (pP38T180/Y182) increased over a broader concentration range, and the weak increase in activation loop–phosphorylated ERK (pERKT202/Y204) slowly approached saturation. The Hill exponent for JNK activation, which quantifies the ultrasensitivity of the response, was greater than that observed for either p38 or ERK activation (Fig. 1C). The abundance of total JNK, p38, and ERK did not change under these conditions (fig. S1D). Additionally, a p38 inhibitor had no effect on apoptosis induced under any of the tested conditions (fig. S1E), indicating that switch-like activation in response to stress was a distinctive aspect of the apoptotic JNK pathway in SH-SY5Y cells under these conditions.

Identifying the components of the JNK network in neuroblastoma cells

The network structure that mediates this behavior of JNK is uncharacterized in neuroblastoma but could arise through positive feedback from JNK to upstream components of the pathway. To investigate the existence of JNK-mediated positive feedback, we experimentally disrupted that loop using a JNK inhibitor and measured the activation loop phosphorylation of MKK4 and MKK7 in response to increasing anisomycin concentration (Fig. 1D). In the presence of a JNK inhibitor, activation of JNK and MKK7 (detected as the forms phosphorylated in the activation loops) by anisomycin was significantly reduced and no longer resembled a high-amplitude ultrasensitive response, whereas the activation loop phosphorylation of MKK4 was less affected. The abundance of MKK4, MKK7, and JNK did not change under these conditions (fig. S1D). We also observed similar changes in response to anisomycin when the cells were exposed to the structurally unrelated JNK inhibitor VIII (fig. S2A) or in cells exposed to the osmotic stress agent sorbitol in the presence or absence of JNK inhibitor II (fig. S2B). This positive influence of JNK activity on its own upstream kinases conforms to our hypothesis of a JNK-mediated positive feedback loop that is activated during cell stress.

To ascertain the contribution of MKK4 and MKK7 (1) to JNK activation in neuroblastoma cells, we individually knocked down MKK4 and MKK7 by small interfering RNA (siRNA) and analyzed anisomycin-induced JNK activation (fig. S2C). Under these conditions, MKK4 contributed ~30% of JNK phosphorylation and MKK7 contributed ~70%, with the latter siRNA providing a similar amount of inhibition as direct JNK inhibition (compare fig. S2C and Fig. 1D, first graph).

To develop a model of the JNK network, we needed to identify the MAPKKK upstream of MKK4 and MKK7. ZAK is an MAPKKK that mediates anisomycin stimulation of the JNK pathway in other cell types (24, 25). ZAK knockdown in SH-SY5Y cells inhibited the activation of JNK, MKK4, and MKK7 by anisomycin and other stressors (Fig. 1E) but did not alter the total amount of JNK, MKK4, or MKK7 (fig. S2D). Stable ZAK knockdown also increased the concentration at which 50% of the cells died [median inhibitory concentration (IC50)]: ~3-fold for vincristine and ~10-fold for doxorubicin (fig. S2E), further suggesting that the ZAK-MKK4/MKK7-JNK signaling axis is central to the cytotoxic response of neuroblastoma cells to therapeutic agents.

Scaffolding proteins, such as the JNK-interacting protein (JIP) family, can contribute to switch-like JNK signaling through the colocalization of MAPKKK-MAPKK-MAPK modules (26). However, we could not detect JIP1 or JIP4 in SH-SY5Y cells, and knocking down JIP2 or JIP3 by siRNA had no effect on anisomycin-induced JNK activity (fig. S3A). Although we did not detect a role for the JIP scaffold family, we investigated protein-protein interactions within this kinase network. ZAK can activate JNK through both MKK4 and MKK7 (27). After activation, dimerization, and autophosphorylation, ZAK interacts directly with MKK7 (28).

To test for protein interactions, we expressed tagged ZAK with either tagged JNK2 or tagged MKK7 in SH-SY5Y cells. Although we coimmunoprecipitated enhanced green fluorescent protein (EGFP)–tagged ZAK and endogenous MKK7 (fig. S3B), MKK4 and EGFP-ZAK did not coimmunoprecipitate (fig. S3B). JNK tagged with the V5 epitope coiummunoprecipitated with EGFP-ZAK (fig. S3C). Knockdown of MKK7 did not reduce the coimmunoprecipitation of EGFP-ZAK and JNK-V5 (fig. S3D). JNK docks to substrates by binding D domains. ZAK has many putative D domains [based on the (K/R)1-3X1-6(L/I)X(L/I) consensus sequence for JNK binding (29)], and one or more of these putative D domains may enable ZAK to interact directly with JNK. Given that we did not detect any involvement of JIPs, the coimmunoprecipitation data suggest the possibility that ZAK may act as a scaffold for MKK7 and JNK.

Identifying positive feedback from JNK to MKK7

Because the known JNK scaffolds were not involved in JNK activation in this system, we tested whether JNK-mediated feedback phosphorylation of its upstream kinases was responsible for ultrasensitive JNK activation. We used MAPK substrate motif antibodies (pTP and PxpSP or pSPxR/K) to investigate potential JNK-dependent phosphorylation events on MKK4, MKK7, and ZAK by Western blotting. Indeed, EGFP-MKK4 and EGFP-MKK7 underwent anisomycin-stimulated and JNK-dependent phosphorylation at TP motifs, and MKK4 was also phosphorylated at PxSP or SPxR/K motifs (Fig. 2A). The amount of ZAK phosphorylated at TP motifs was unaffected by anisomycin. Because the observed positive feedback loop appeared to act predominantly on MKK7 (Fig. 1D), we focused on feedback phosphorylation of MKK7 by JNK. We mutated Thr66 and Thr83 within TP motifs in MKK7 because they flank a D domain, through which JNK docks to phosphorylate substrates (29). Both mutations significantly impaired anisomycin-induced phosphorylation of MKK7 at TP motifs, with the T66A mutation having a larger effect (Fig. 2B).

Fig. 2 Positive feedback in the JNK network.

(A) Western blotting of immunoprecipitated (IP), GFP-tagged ZAK, MKK4, or MKK7 transfected into SH-SY5Y cells, which were subsequently treated with anisomycin (300 nM, 30 min) in the presence or absence of JNK inhibitor II. Blotting was performed using MAPK substrate antibodies directed toward the PxSP_SPxR/K or TP motifs or antibodies recognizing GFP. (B) Western blotting of immunoprecipitated, GFP-tagged wild-type (WT) MKK7, MKK7T66A, or MKK7T83A transfected into SH-SY5Y cells that were treated with anisomycin (300 nM, 30 min). (C) TP motif phosphorylation of immunoprecipitated, FLAG-tagged murine MKK7β1 and MKK7α1 treated with anisomycin (300 nM, 30 min) in the presence and absence of SP600125 (JNK inhibitor II). IgG, immunoglobulin G. (D) Activation loop phosphorylation of immunoprecipitated, V5-tagged WT MKK7, MKK7T66A, or MKK7T83A transfected into SH-SY5Y cells that were treated with anisomycin (300 nM, 30 min). In panels (A), (B), and (D), quantification is from three independent experiments (mean ± SD, n = 3). *P < 0.05; **P < 0.01, two-sample Student’s t test. Phosphorylation is calculated relative to the untreated control sample, which was set at 1. (E) Analysis of the coimmunoprecipitation of GFP-tagged ZAK and V5-tagged MKK7 (WT, T66A, and T83A) after anisomycin treatment in SH-SY5Y cells. Data are representative of three experiments.

To further investigate the mechanism underlying this JNK-mediated positive feedback, we used the murine MKK7β1 isoform, in which Thr66 and Thr83 are conserved, and the murine MKK7α1 isoform, which lacks the N-terminal region (Ala2-His88) that contains these feedback phosphosites. We confirmed that FLAG-tagged murine MKK7β1, but not MKK7α1, was phosphorylated in a JNK-dependent manner in response to anisomycin when expressed in SH-SY5Y cells (Fig. 2C). Mutation of Thr66, but not Thr83, within MKK7β1 significantly reduced anisomycin-induced MKK7 activation loop phosphorylation (Ser271), confirming that this site constitutes the positive feedback from JNK to MKK7 (Fig. 2D). In contrast, the MKK7α1 isoform, which cannot undergo JNK-mediated feedback phosphorylation, displayed significantly increased basal activation loop phosphorylation, which did not change upon anisomycin treatment (Fig. 2D).

Another possible mechanism for the feedback could be an enhanced interaction between ZAK and MKK7. However, the coimmunoprecipitation between ZAK and MKK7 was retained in the MKK T66A or T83A mutants, indicating that feedback phosphorylation from JNK to MKK7 did not affect the binding of MKK7 to ZAK (Fig. 2E). These protein interaction data and the phosphorylation analysis indicated that the N-terminal region of MKK7 may be an autoinhibitory domain and that inhibition is relieved upon JNK-mediated phosphorylation, thereby providing a mechanistic basis for the observed positive feedback loop within the JNK network.

AKT-mediated crosstalk to JNK

Another network connection capable of regulating JNK activation is AKT-mediated inhibitory phosphorylation of MKK4 at Ser80 (30). To explore the crosstalk between AKT and JNK signaling in SH-SY5Y cells, we transfected a constitutively active, myristoylated AKT1 (Myr-AKT1) and examined the phosphorylation of the kinases in the JNK pathway. Myr-AKT1 increased phosphorylation of MKK4 at Ser80, correlating with the reduced activation loop phosphorylation of MKK4 and JNK in response to anisomycin and vincristine (Fig. 3A). Moreover, Myr-AKT1 also decreased the stress-induced activation loop phosphorylation of MKK7.

Fig. 3 AKT inhibits JNK activation.

(A) Western blotting of lysates from SH-SY5Y cells transfected with Myr-AKT1 or a control plasmid and treated with anisomycin (300 nM, 30 min) or vincristine (100 nM, 1 hour). (B) Western blotting of immunoprecipitated, GFP-tagged ZAK, GFP-tagged MKK4, or GFP-tagged MKK7 cotransfected into SH-SY5Y cells with Myr-Akt1 or a control plasmid. Blotting was performed using antibodies specific for Akt substrate phosphorylation (RxRxxT/S_RxxT/S) or GFP. (C) Western blotting of lysates and immunoprecipitated, GFP-tagged MKK7 or GFP-tagged MKK7T385A cotransfected into SH-SY5Y cells with Myr-Akt1 or a control plasmid. Quantification is from three independent experiments (mean ± SD, n = 3). *P < 0.05; **P < 0.01, two-sample Student’s t test. Phosphorylation was calculated relative to the untreated control sample, which was set at 1.

To detect other AKT substrates in this MAPK cascade, we cotransfected Myr-Akt1 and GFP-tagged MKK4, GFP-tagged MKK7, or GFP-tagged ZAK; immunoprecipitated the GFP-tagged proteins; and assayed the proteins for phosphorylation with an antibody that recognizes the phosphorylated AKT consensus site (RXRXXpS/pT or RXXpS/pT) (Fig. 3B). Both MKK4 and MKK7 exhibited increased phosphorylation in cells cotransfected with Myr-AKT1 compared with the amount in cells expressing only the GFP-tagged kinase. However, phosphorylation of ZAK was the same in cells cotransfected with GFP-ZAK and Myr-AKT1 and in cells expressing only GFP-ZAK, indicating that phosphorylation recognized by this antibody was not mediated by AKT. MKK7 does not contain any canonical RXRXXS/T AKT phosphorylation motifs; however, it does contain an RXXT motif (381RYET385), which is sufficient for Akt-mediated phosphorylation of threonine (31). Mutation of Thr385 to alanine resulted in the loss of AKT-mediated phosphorylation of MKK7 (Fig. 3C), confirming that both MKK4 and MKK7 are AKT substrates.

Mathematical modeling of the JNK network dynamics

To generate quantitative predictions using this experimentally resolved network structure, we built, calibrated, and validated a computational model of JNK network dynamics (Fig. 4A). This model is based on ordinary differential equations (ODEs) and incorporates the following experimentally confirmed features. First, stress activates ZAK, which activates MKK4 and MKK7, which in turn activate JNK. Second, there is a positive feedback from JNK to MKK7 and an inhibitory crosstalk from AKT to MKK4 and MKK7. Third, ZAK functions as both a kinase and a scaffold that binds MKK7 and JNK. The scaffold-mediated interactions, double reversible phosphorylation cycles in the MAPK cascade, and AKT crosstalk result in a combinatorial increase of feasible reactions and multiple states of the JNK network kinases (32, 33). Therefore, we built the model using rule-based modeling (3437), which develops the ODEs on the basis of a framework of rules describing the reactions and states of the system. Constructed in this fashion, the model implemented 23 rules, which resulted in 258 reactions and 36 parameters (table S1). The rule-based approach minimizes the number of necessary simplifying assumptions, and although it increases the number of states and reactions, it does not increase the number of parameters in the model.

Fig. 4 Modeling of the JNK network dynamics.

(A) Graphical representation of the model as a limited process graph (for clarity, multiprotein complexes, double-phosphorylation cycles, and multiple inhibitory forms are absent). Asterisks (**) indicate double-phosphorylated, active forms of the kinases. (B) Model fitting to experimental data. Lines show the average of model simulations for all 50 parameter estimates (mean ± SD, bold line and shaded areas, respectively) and the best-fitting model (dashed lines). (C) Comparison of simulated (left) and observed (middle and right) anisomycin-induced pJNKT183/Y185 in SH-SY5Y cells after knockdown of MKK4 or MKK7 and overexpression of ZAK, as indicated. The arrow indicates the band specific for MKK4. The three bands for JNK represent different splicing variants of JNK1, JNK2, and JNK3; the intensities of all three bands were included for quantification. (D) Prediction of the stress-induced phosphorylation of JNK in the indicated neuroblastoma cell lines. Left panel: Representative blots and measurements (mean, n = 3 blots) of relative protein abundance, normalized to actin. Unless marked by an arrow, all bands were used for quantification. Right panel: Simulated activation of JNK based on these parameters. (E) Experimental validation of the predicted JNK responses in panel (D). Left panel: Activation loop phosphorylation of JNK in the three cell lines after vincristine treatment (100 nM, 2 hours). Right panel: Linear correlation between predicted and observed pJNKT183/Y185, normalized to the loading control. The small red dots represent the best-fitting model. (F) Linear correlation (dashed line) between experimental (observed) or simulated (predicted) pJNKT183/Y185 and apoptosis in response to vincristine (100 nM, 24 hours). In panels (E) and (F), experimental data are mean ± SD, n = 3; simulations are mean ± SD, n = 50. siContr., control siRNA; siMKK4 or siMKK7, siRNA for MKK4 or MKK7; Ani, anisomycin; Unstim., unstimulated (not exposed to stress-inducing drug).

To generate quantitative predictions, we calibrated the model using our time course (fig. S1C), dose-response (Fig.1C), and MAPKK knockdown data (Fig. 1E). We adjusted key parameters, such as the kinase activities and feedback strengths, using a global parameter estimation method called adaptive simulated annealing (38). To assess the associated uncertainty, we adopted a Monte Carlo–based approach by randomly changing the initial parameters 50 times and refitting the model (table S2 and text S1) (39, 40). We then analyzed parameter correlations and variability of the simulated predictions. Although parameter estimates have large SDs, the estimates do not spread over the entire multidimensional space but occupy structured regions where some parameter values are correlated (fig. S4, A to C). Although exact parameter values cannot be identified, the predictions are accurate and tight for all these estimates (Fig. 4C).

Experimental model validation

Model simulations performed to simulate the presence of a JNK inhibitor recapitulated our experimental observations of a positive feedback loop, in that blocking the positive feedback from JNK to MKK7 impaired the stress response and its switch-like, ultrasensitive behavior (Fig. 4B, middle column). We further tested the predictive power of the model by examining the potential scaffolding effect of ZAK (Fig. 4C, left), using an independent validation data set (Fig. 4C, middle and right). Simulations predicted that ZAK overexpression would impair JNK activation because of the inability of ZAK, when in excess, to simultaneously scaffold both MKK7 and JNK. Additionally, these simulations predicted that ZAK overexpression would further impair JNK activation after MKK4, but not MKK7, knockdown. In support of our experimentally resolved network structure, these model predictions were fully corroborated by experimental analysis of anisomycin-mediated JNK activation after the simultaneous overexpression of ZAK and knockdown of either MKK4 or MKK7 (Fig. 4C).

Initially, our model calibration and validation experiments were performed in SH-SY5Y cells, using anisomycin as a stressor. We therefore tested whether the model correctly predicted the JNK response to different stress-inducing agents and in other neuroblastoma cell lines. To this end, we measured the abundance of each kinase included in the model and the phosphorylation of AktS473 in the SMS-KCN and IMR32 lines relative to their amounts in the SH-SY5Y line (Fig. 4D, left). We then populated our original model with these values as parameters, thereby generating three cell line–specific models. Simulations predicted that the JNK response would be markedly impaired in SMS-KCN cells and almost absent in IMR32 cells (Fig. 4D, right). The predicted differences in the JNK pathway activity (JNK activation loop phosphorylation) among the three cell lines in response to stress were confirmed experimentally for treatment with vincristine (Fig. 4E), thus demonstrating that our model can be used to simulate and correctly predict the JNK response in different neuroblastoma cell lines and across different stress stimuli. Furthermore, both the observed and predicted JNK responses for each cell line were proportional to the apoptotic response induced by vincristine (Fig. 4F), suggesting that our model also has predictive power with regard to the biological response of neuroblastoma cells to clinically relevant therapeutic agents.

Prognostic utility of the JNK model

On the basis of our combined experimental and modeling studies, we hypothesized that a high-amplitude ultrasensitive JNK stress response would promote apoptosis in neuroblastoma cells and would therefore be associated with a better patient prognosis than a low-amplitude or gradual JNK response. To test this hypothesis, we generated patient-specific models by incorporating gene expression data from three cohorts of primary neuroblastoma samples (table S3): a training cohort of 109 patients and two independent validation cohorts of 369 and 233 patients. Although the correlation between mRNA and protein abundance is weakly positive on a global scale (41), our analysis of multiple data sets consistently provided a strong correlation for all proteins in our model (fig. S5, A to C). Indeed, each of the transcripts and proteins exhibited correlation coefficients ≈0.6 and proportionality factors ≈1 (text S2 and tables S4 and S5). Therefore, we adjusted the protein concentrations in the model for each patient according to the relative gene expression measured in the tumor (fig. S6A, Fig. 5A, and fig. S7A); for example, if the tumor sample showed a twofold increase in MKK4 mRNA, then we increased the abundance of MKK4 by twofold in the model. Following this procedure for all model components and for each patient resulted in personalized models that incorporated protein abundance on the basis of the measured gene expression of ZAK, MKK4, MKK7, JNK (average of MAPK8, MAPK9, and MAPK10 probes), and AKT (average of AKT1 and AKT2 probes).

Fig. 5 Patient-specific simulations and model-based stratification.

(A) Heat map representing the relative gene expression (mRNA) within a cohort of 369 neuroblastoma patients. Each column corresponds to one patient. (B) Simulated JNK stress response for each individual patient using the data in (A) as parameters in the model. Thin lines are the individual responses; bold lines are the means for each stage. (C) Network output descriptors derived from the patient-specific simulations: amplitude (A), half-activation threshold (K50), and Hill exponent (H). (D) Distribution of the simulated network output descriptors for the different stages. (E) Model-based Kaplan-Meier survival analysis for the entire cohort of 369 patients. Top panel: Scanning for the best separation cutoff; the red dot indicates the optimized cutoff derived from the calibration cohort (fig. S6). The statistical significance (log-rank test, y axis) of the overall survival difference is plotted as a function of the cutoff [number of patients in the low group (%), x axis]. Middle panel: The population’s network output descriptors profile illustrating the optimal low and high groups. Values represent the optimal cutoff values for A, K50, and H. Bottom panel: Kaplan-Meier overall survival curves for the low and high groups (P value log-rank test, with Yates’ correction). HR, hazard ratio (log2 scale); +, censored data for patients without death event; iqr, interquartile range; w.r.t., with regard to.

Having generated these personalized models of the JNK stress response for each patient, we characterized these dynamic responses by extracting three JNK network output descriptors: the maximal amplitude (A), the activation threshold (K50), and the Hill exponent (H) (fig. S6, B and C, Fig. 5, B and C, and fig. S7, B and C). Using these output descriptors, an impaired ability to activate JNK corresponded to high values for activation threshold and low values for maximal amplitude and Hill exponent. Grouping these patient-specific simulations according to the International Neuroblastoma Staging System (INSS) revealed increasingly impaired JNK responses as the disease progressed (fig. S6, B and D, Fig. 5, B and D, and fig. S7, B and D). This impairment was apparent from the decreased mean response in higher disease stages, described by significantly decreased amplitude in stages 4 and 4S, significantly increased activation threshold in stages 3 and 4, and a significantly reduced Hill exponent in stage 4 and 4S.

To further examine the association between our model of JNK activation and disease progression, we stratified the patient population into two groups according to their simulated JNK network output descriptors. To identify the optimal cutoff values that best correlate or anticorrelate with survival for each output descriptor, we applied Kaplan-Meier scanning (42) to the training cohort (fig. S6E). Briefly, we sorted the patients according to each output descriptor and started by putting 10% of patients into the low group and 90% into the high group and then calculating the log-rank test P value in iterative steps. In each step, the size of the low group is increased by one patient. The best cutoff is the value of the output descriptor with the lowest P value. This patient stratification process establishes the values that can be used to statistically significantly associate survival with Hill exponent, maximal amplitude, or threshold for JNK activation. Plotting the corresponding Kaplan-Meier survival curves showed a significant association between poor patient survival and low amplitude, low Hill exponent, or high threshold of JNK activation (fig. S6E). We then applied the same cutoffs to the independent validation cohorts and found the same association between an impaired JNK response and poor overall survival (Fig. 5E and fig. S7E). Notably, the Hill exponent, which measures the ultrasensitivity of the JNK response, provided the most accurate patient stratification. These results suggest that impairment of the JNK stress response, resulting from alterations in the concentration of kinases within the JNK network, contributes to a negative outcome for neuroblastoma patients.

JNK network alterations induced upon MYCN amplification

Because MYCN is the major oncogene driving neuroblastoma (19), we hypothesized that MYCN-amplified tumors may have altered JNK stress response dynamics. Inspection of each kinase in our model revealed that the most prominent changes in MYCN-amplified tumors were decreased MKK4 and increased ZAK and AKT (Fig. 6A), all of which would be expected to impair JNK responses according to our model and experimental observations. Indeed, model simulations showed that MYCN-amplified patients exhibited markedly impaired JNK responses when compared to patients without MYCN amplification, with significant alterations in all three output descriptors of the JNK network (Fig. 6, B and C). These coordinated changes in the transcripts for the JNK network kinases and their association with MYCN amplification suggest that an impaired JNK response may be characteristic of MYCN-amplified tumors and may contribute to the innate chemoresistance associated with this aggressive form of neuroblastoma (19, 43).

Fig. 6 MYCN-associated changes and MYCN-independent effects.

(A) Distribution of measured mRNA abundance for JNK network components within the non–MYCN-amplified (0) and MYCN-amplified (1) cohorts (**P < 0.001; *P < 0.05, two-sample Student’s t test). (B) Individual simulated stress responses (thin lines) for both cohorts and the corresponding mean responses (thick black lines). (C) Box plots and group differences for the simulated network output descriptors calculated from (B). (D) Model-based Kaplan-Meier survival analysis for the non–MYCN-amplified and MYCN-amplified cohort. See Fig. 5E for a description of all elements.

Despite this potential for MYCN to drive changes that impair the JNK response, a substantial number of patients lacking MYCN amplification also exhibited an impaired JNK response in our simulations. We therefore asked whether our model could also predict possible MYCN-independent effects on overall survival by performing survival analysis for the non–MYCN-amplified and MYCN-amplified groups separately. Under these conditions and considering P < 0.01 significant, only the Hill exponent was significantly associated with overall survival for both the MYCN-amplified and non–MYCN-amplified patients when using cutoffs derived from the calibration cohort (Fig. 6D and fig. S8, A and B). However, when the optimal cutoffs were reevaluated for each subcohort (MYCN-amplified and non–MYCN-amplified) separately, all calculated network output characteristics of an impaired JNK response were significantly associated with poor patient survival in the nonamplified cohort (fig. S8A). For MYCN-amplified patients, the Hill exponent was significantly associated with poor survival (fig. S8B), and the amplitude (P = 0.026) came close to reaching our significance threshold (P < 0.01). Thus, our modeling approach has confirmed that an impaired JNK response is significantly associated with poor survival in both MYCN-amplified and non–MYCN-amplified cases.

We further tested whether using the model improves the current prognostic understanding based on well-established markers, such as MYCN status, age, and stage (INNS) or risk [Children’s Oncology Group (COG)]. Using several Cox proportional hazard analyses correcting for each marker in different combinations, we found that the dynamic model reliably contributed to predicting the overall survival (Table 1). Individually as independent variables in a univariate analysis, maximal amplitude, Hill exponent, and threshold for JNK activation all significantly predicted overall survival. When combined with other markers in multivariate analysis, only maximal amplitude and Hill exponent, but not activation threshold, were independent prognostic indicators.

Table 1 Cox proportional hazard analysis of overall survival using the JNK output descriptors as key independent variables.

Univariate analysis was performed using a Cox proportional hazard model with either A, K50, or H as an independent variable. Multivariate analysis was performed using a Cox proportional hazard model with the indicated variables measured at the time of admission or surgery as covariates in addition to either A, K50, or H. β is presented as the maximum likelihood estimate ± SE. Risk is defined by COG. MYCN indicates positive or negative for MYCN amplification.

View this table:

Robustness of the JNK dynamics and prognostic predictions

To determine whether the network model increased the amount of obtainable information or merely identified useful individual biomarkers, we compared the prognostic utility of the individual model components to that of the network-based model. Performing survival analyses on each model component in isolation revealed significant associations between overall survival and the expression of ZAK, MKK4, AKT (P < 0.00001), and JNK (P < 0.01), but not MKK7 (P ≈ 0.05) (fig. S9A). In the cohort without MYCN amplification, increased ZAK and decreased MKK4 were significantly associated with poor overall survival (P < 0.01) (fig. S9B). In the MYCN-amplified cohort, no single model component was significantly associated with survival (fig. S9C). MKK7 (P = 0.02443) replaced MKK4 (P = 0.03742) as the most useful individual prognostic indicator in MYCN-amplified patients.

MKK4 showed the strongest correlation with overall survival and was more significantly associated with poor overall survival than our modeling analysis in the whole patient cohort. We therefore investigated the contribution of MKK4 to our modeling analysis in more detail. A sensitivity analysis of the model confirmed that changing the concentration of MKK4 markedly affected the JNK response (Fig. 7A). Although MKK7 and JNK were the most sensitive parameters with respect to equally sized perturbations (Fig. 7A, left, and fig. S10A), MKK4 was the most sensitive when the perturbation size was corrected for the variability observed in the tumor data (Fig. 7A, right) and had the predominant effect on shaping both the JNK time course and dose response (fig. S10B). This finding suggests that although MKK4 does not play a dominant role within this network structure, the abundance of MKK4 in the tumors is more markedly altered than that of the other model components. Accordingly, MKK4 contributed independently in a multivariate Cox regression analysis using MYCN amplification, age, and all kinases in the model as covariates (Table 2). However, whereas MKK4 expression was associated with survival in the whole patient cohort, when analyzing only MYCN-amplified or high-risk patients, the model outperformed the predictive capabilities of MKK4 (Fig. 6B). Indeed, in the MYCN-amplified cohort, the model outperformed the predictive capacity of any of the individual components in the model (compare Fig. 5D and fig. S9C).

Fig. 7 Robustness of the JNK dynamics.

(A) Sensitivity analysis and perturbation analysis of the JNK network. (Left) The effect on JNK activation of a onefold perturbation in the concentration of each kinase in the model. (Right) The effect of perturbing the kinase concentration according to the measured variation within the tumor data. For each kinase, the SD in the tumor data was calculated, and the amount of the corresponding component in the model was increased (Up) or decreased (Down) by twice that value. Perturbation values are indicated in parentheses after the name of each kinase. (B) Receiver operating characteristic (ROC) analysis to determine the predictive quality of the model, MYCN status, or abundance of MKK4. The patient outcome (death or survival) was predicted using the simulated H from the model or the measured MYCN mRNA or MKK4 mRNA expression. The resulting true-positive rate (y axis) was plotted against the false-positive rate (x axis). AUC denotes the area under the curve. (C) Distribution of the network output descriptors across the whole cohort (n = 369) with and without positive feedback from JNK to MKK7. (D) Distribution of the network output descriptors across the disease stages, performed without positive feedback in the model (**P < 0.001; *P < 0.05, two-sample Student’s t test). Pos. FB, positive feedback; No FB, no feedback; CV, coefficient of variation defined as the SD divided by the mean.

Table 2 Multivariate Cox regression analysis.

A Cox proportional hazard model was used to analyze the overall survival for the entire cohort of all patients.

View this table:

These cohort-specific associations suggest that placing these individual signaling pathway components into their network context is a necessary step for obtaining personalized simulations that can deliver robust predictions against different molecular backgrounds. This importance of network wiring was particularly highlighted by performing the patient simulations after removal of the positive feedback from within our model (fig. S10C). This analysis resulted in a substantial decrease in the heterogeneity of the simulated patient responses for both the amplitude and Hill exponent (Fig. 7C) and a complete loss of the prognostically relevant variability across the disease stages for the K50 values and the Hill exponent (Fig. 7D). The loss of stage-specific variation in JNK ultrasensitivity is of particular note given that the Hill exponent provided the most accurate stratification of the MYCN-amplified cohort. Further, it suggested that the JNK-MKK7 positive feedback enhances the overall impact of small alterations on the abundance of network components, thereby boosting the dynamic network response and leading to the generation of biologically meaningful response changes that ultimately influence patient survival.


Our experimental characterization and personalized modeling of stress-induced JNK activation highlights the central role of JNK signaling in neuroblastoma patient outcome. The identification of major genetic alterations in JNK network components in breast and pancreatic cancer (4446) further suggests that preventing JNK activation is a vital step during tumor progression and may ultimately be achieved through diverse mechanisms. Here, we identified mechanisms through which alterations in the transcript abundance and, thus, protein abundance of JNK network components result in impaired ultrasensitive JNK activation in neuroblastoma cells. The importance of these mechanisms is also supported by other studies that identified the ZAK-MKK4/MKK7-JNK module as a crucial mediator of apoptotic signaling in multiple cell types (24, 25, 47). Indeed, studies have implicated the off-target inhibition of ZAK by inhibitors of the MAPKKK BRAF as a mechanism of chemoresistance in BRAF-mutant melanoma cells (48, 49).

The importance of adopting a modeling approach when investigating the connections within a signaling network and the implications for disease is highlighted by the prognostic significance of the patient-specific simulations from tumors with MYCN amplification. In this cohort, the simulation of JNK activation was significantly associated with overall survival, whereas none of the survival analyses performed with individual model components reached significance. Further refinement of this patient-specific modeling approach could certainly be achieved by using protein-based data [for example, reverse-phase protein array (RPPA)], particularly regarding the inclusion of direct measurements of AKT phosphorylation. However, the ability of the model to achieve prognostic significance based on simulations performed using gene expression data lends credence to the strength of this modeling approach and is also more practical for clinical application because mRNA expression analysis is more common than RPPA.

The identification of a positive feedback loop from JNK to MKK7, facilitating the switch-like activation of JNK in our model, provides a potential mechanism for many experimental observations of ultrasensitive JNK activation (12, 14, 15). Additionally, the combination of this positive feedback and patient-to-patient variation in the expression of transcripts encoding the JNK network components also underpins the prognostically significant heterogeneity in the patient-specific simulations JNK activation. So far, it has been unclear how quantitative changes in expression might contribute to response variability, cell fate, and survival because small changes might be statistically significant but not biologically meaningful. Here, we demonstrated that many small changes distributed across the kinases within the JNK pathway synergize to evoke biologically meaningful response differences that translate into differences in patient survival. Small variations at the level of gene expression are enhanced at the level of pathway activation by the positive feedback that we identified within the JNK pathway. The response heterogeneity introduced by the positive feedback from JNK to MKK7 contrasts with the related ERK pathway, in which negative feedback from ERK to RAF-1 increases robustness to perturbations and fluctuations in protein abundance (50, 51). These findings highlight the vital role of compiling detailed information about network wiring in the development of computational models of signaling pathways. Even highly related pathways are not necessarily wired in a similar fashion, and any assumptions need to be carefully validated to derive biochemically correct and clinically relevant information.

Computational modeling has already provided an understanding of the dynamics of signaling pathways (16, 5254). Here, we showed that modeling approaches can be used not only to decipher and understand the behavior of complex networks on the cellular level but also to relate biochemical pathway behavior to clinical outcomes on the level of individual patients. The simulation of patient-specific apoptotic JNK network responses provides an important proof of principle for the potential prognostic utility of experimentally resolved models of signaling networks, as opposed to biomarkers that lack this detailed mechanistic anchoring. The crucial advantage of a mechanistic model in this setting is that it can be quickly adapted to compare and predict the efficacy of different therapeutic agents once their mode of action has been determined.


Study design

This study used the SH-SY5Y, SMS-KCN, and IMR32 neuroblastoma cells lines to generate, calibrate, and validate a mathematical model of stress-induced JNK signaling. Patient-specific simulations of JNK signaling were performed using three existing independent patient cohorts (55). The training cohort included 109 international patients, the first validation cohort included 369 German patients, and the second validation cohort included 233 international patients (table S3).

For the training cohort, patients’ age at diagnosis ranged from 3 days to 293 months (median age, 12.7 months), median follow-up for patients without fatal events was 5.4 years (range, 0.2 to 18 years), and 5-year overall survival was 0.63 (95% CI, 0.53 to 0.73). For the German validation cohort, patients’ age at diagnosis ranged from 0 to 295 months (median age, 13.2 months), median follow-up was 3.4 years (range, 0 to 16 years), and 5-year overall survival was 0.80 (95% CI, 0.75 to 0.86). For the second validation cohort, patients’ age at diagnosis ranged from birth to 304 months (median age, 17 months), median follow-up for patients without fatal events was 4.5 years (range, 0.02 to 24 years), and 5-year overall survival was 0.73 (95% CI, 0.67 to 0.80). Stage was classified according to the INSS, and risk according to the COG system.

Antibodies, plasmids, and stress-inducing reagents

Antibodies specific to total JNK (#9258), MKK4 (#9152), MKK7 (#4172), pJNKT183/Y185 (#9251), pERKT202/Y204 (#4370), pP38T180/Y182 (#4511), pMKK4S257/T261 (#9156), pMKK4S80 (#9155), pAktS473 (#4060), pMAPK substrate antibodies T*P (#9391) and PxS*P_S*PxR/K (#2325), pAkt substrate RxRxxT/s, RxxT/S (#9614), and cleaved caspase 3 (#9664) were from Cell Signaling. The pMKK7S271 (TA312581) antibody was from OriGene. The GFP monoclonal antibody (11814460001) was from Roche. The antibodies toward ZAK (sc134970), p53 (sc126), JIP1 (sc25267), JIP2 (sc53553), and JIP4 (sc271491) were from Santa Cruz Biotechnology. The JIP3 antibody (NBP1-00895) was from Novus Biologicals. The actin monoclonal antibody (AC-15) and FLAG-M2 antibody were from Sigma-Aldrich. The GAPDH (ab8245) and MKK4 (ab33912) antibodies were from Abcam. The V5 monoclonal antibody (R960-25) was from Life Technologies.

The pDONR223-MKK7 construct (56) was a gift from W. Hahn and D. Root (Addgene plasmid 23798). The pDONR201-MKK4 construct (HsCD00000136) (57) was obtained through DNASU (58). The pDONR-ZAK construct (GC-X0399-CF) was purchased from GeneCopoeia. Cloning of these donor cDNAs into 223 pCS EGFP and pcDNA6.2 V5 destination plasmids was performed using the Gateway LR Clonase enzyme mix (Life Technologies) according to the manufacturer’s instructions.

JNK inhibitor II (SP600125) was from Merck, and JNK inhibitor VIII was from Cayman Chemical. Anisomycin, sorbitol, rotenone, vincristine, and doxorubicin were from Sigma-Aldrich. Anisomycin dose responses were performed using 3, 10, 30, 100, 300, and 1000 nM anisomycin for 30 min. Sorbitol dose responses were performed using 10, 30, 100, 300, 500, and 1000 nM sorbitol for 1 hour. Rotenone dose responses were performed using 3, 10, 30, 100, 300, and 1000 nM rotenone for 1 hour. Vincristine dose responses were performed using 3, 10, 30, 100, 300, and 1000 nM for 2 hours. Doxorubicin dose responses were performed using 1, 3, 10, 30, 100, and 300 nM for 4 hours.

Zebrafish model of neuroblastoma

The transgenic Tg[dβh:MYCN-EGFP] zebrafish line was a gift from A. T. Look (Dana-Farber Cancer Institute, Harvard Medical School) (23). Zebrafish experiments were approved by the University College Dublin Animal Research Ethics Committee and the Healthcare Products Regulatory Authority (AE18982/P038).

Adult Tg[dβh:MYCN-EGFP] zebrafish (Danio rerio) were maintained at 28°C on a 14-hour light/10-hour dark cycle and monitored for fluorescent tumor masses every 2 weeks. About 20% of fish develop tumors by ~8 months (23). Zebrafish developing the fish-equivalent of human neuroblastoma (age, 15 months; n = 4) were treated every 48 hours with JNK inhibitor II (5 μM) or DMSO control for a total of 2 weeks. Variance in tumor size for all fish at pretreatment was <5%. Zebrafish were anesthetized by 0.016% tricaine before microscopic analysis. An Olympus SZX10 fluorescence stereomicroscope equipped with an Olympus DP71 camera was used to capture images. The fluorescent region posterior to the gills, which corresponds to the region of neuroblastoma growth (23), was analyzed using ImageJ software (version 1.44p).

Cell lines

The SH-SY5Y, SMS-KCN, and IMR-32 neuroblastoma cell lines have been described previously (59). These cell lines were cultured in standard RPMI 1640 (containing d-glucose at 2 g/liter) under standard tissue culture conditions (5% CO2, 20% O2).

Site-directed mutagenesis

Site-directed mutagenesis of MKK7 was performed using the QuikChange Lightning Site-Directed Mutagenesis Kit (Agilent) according to the manufacturer’s instructions. The following primers were used for human MKK7: T66A forward primer, 5′-gcagcaccccgcgccccccgc-3′, and T66A reverse primer, 5′-gcggggggcgcggggtgctgc-3′; T83A forward primer, 5′-catgctgcggggtgcgaacagggttgacg-3′, and T83A reverse primer, 5′-cgtcaaccctgttcgcaccccgcagcatg-3′; T385A forward primer, 5′-caagcgctacgaggcgctggaggtgga-3′, and T385A reverse primer, 5′-tccacctccagcgcctcgtagcgcttg-3′. The following primers were used for murine MKK7: T66A forward primer, 5′-ccacagcaccctgcaccccccaccc-3′, and T66A reverse primer, 5′-gggtggggggtgcagggtgctgtgg-3′; T83A forward primer, 5′-gctcccatcaaccttgttcgcaccgcgcagtat-3′, and T83A reverse primer, 5′-atactgcgcggtgcgaacaaggttgatgggagc-3′.

siRNA and shRNA

The MKK4 (L-003574-00-0005), ZAK (L-005068-00-0005), JIP1 (L-003595-00-0005), JIP2 (L-012462-00-0005), JIP3 (L-003596-00-0005), JIP4 (L-017462-00-0005), and negative control (D-001810-10-05) ON-TARGETplus SMARTpool siRNAs were from Thermo Fisher. The MKK7 Silencer Select siRNA (s11184) was from Ambion. Transfection of siRNA was performed using the jetPRIME reagent (Polyplus Transfection) according to the manufacturer’s instructions. The pLKO control short hairpin RNA (shRNA) and pLKO ZAK shRNA (clone name: NM_016653.x-362s1c1) plasmids were gifts from W. Gallagher (Conway Institute, University College Dublin). Plasmid transfection was performed using the jetPRIME reagent (Polyplus Transfection) according to the manufacturer’s instructions, and stable cell line selection was performed with puromycin (2 μg/ml).

Western blotting and immunoprecipitation

Lysates for Western blotting and immunoprecipitation were prepared using normal lysis buffer [10 mM tris-HCl (pH 7.5), 140 mM NaCl, and 1% Triton X] containing protease inhibitor cocktail (p8340, Sigma) and PhosSTOP Phosphatase Inhibitor Cocktail (Roche). Immunoprecipitation was performed using either GFP-Trap_A (Chromotek), V5 monoclonal antibody–coupled agarose beads, or FLAG-M2–coupled magnetic beads (Sigma) as indicated. SDS–polyacrylamide gel electrophoresis (SDS-PAGE) electrophoresis and Western blotting were performed using the NuPAGE SDS-PAGE gel system and NuPAGE Bis-Tris Precast Gels (10% and 4 to 12%) (Life Technologies). SuperSignal West Femto Chemiluminescent Substrate (Thermo Scientific) was used to develop Western blots, which were imaged using an Advanced Molecular Vision chemiluminescence imaging system. Quantitative Western blotting was performed using multistrip Western blotting (60).

Cell-based assays

Apoptosis was measured by propidium iodide staining and flow cytometric analysis to identify the sub-G1 population, as previously described (61). Analysis was performed using an Accuri C6 flow cytometer (BD). Cytotoxicity assays were performed using the CellTiter 96 AQueous One Solution Cell Proliferation Assay according to the manufacturer’s instructions (Promega).

Dynamic modeling

On the basis of the experimentally resolved network structure depicted in Fig. 4A, a dynamic model was constructed by rule-based modeling (34), which develops the ODEs based on a framework of rules describing the reactions and states of the system. The modeling is described in detail in text S1. An extended contact map provides a graphical representation of this rule-based model (fig. S11). A full account of all rules, equations, and parameters is provided in table S1. The model was implemented in MATLAB using the PottersWheels toolbox, which was also used for parameter estimation using adaptive simulated annealing as an optimizer. Table S2 contains the parameter estimation results for all 50 independent parameter estimation runs. The effect of violating critical model assumptions was also analyzed (fig. S12).

Measuring mRNA expression in neuroblastoma tumor samples

Sample set composition, sample preparation, and generation of single-color gene expression profiles from primary neuroblastomas were described previously (62).

Kaplan-Meier survival analysis

The best cutoff for stratifying the patient population into low and high groups was identified in the training cohort (109 patients) by scanning the group sizes from 10-90 to 90-10 percent splits, where 10-90 means that 10% of the patients were in the low group and 90% of the patients in the high group, and calculating the P value for the overall survival difference between the groups using a log-rank test with Yates’ correction. These optimized cutoff values were then applied to the independent validation cohort (369 patients). Additionally, the optimal cutoffs were also re-derived in the validation set for the separate non–MYCN-amplified and the MYCN-amplified analysis. All statistical computations and Kaplan-Meier analyses were performed in MATLAB using the statistics toolbox and the log-rank ( and kmplot ( functions from the MATLAB file exchange. Multivariate survival analysis is described in text S3.


Text S1. Dynamic modeling of the JNK network.

Text S2. Rationale for using measured mRNA abundances as proxy for protein abundances.

Text S3. Multivariate survival analysis.

Fig. S1. Stress-mediated apoptotic JNK activation.

Fig. S2. Network structures mediating stress-induced JNK activation.

Fig. S3. Detailed wiring of the JNK network.

Fig. S4. Parameter estimation results.

Fig. S5. Correlation between mRNA and protein measurements.

Fig. S6. Patient-specific simulations and model-based stratification in the training cohort.

Fig. S7. Patient-specific simulations and model-based stratification in the second validation cohort.

Fig. S8. Model-based survival analysis for the non–MYCN-amplified and MYCN-amplified cohorts.

Fig. S9. Survival analysis performed using individual model components.

Fig. S10. Robustness of the JNK dynamics.

Fig. S11. Graphical representation of the rule-based model in terms of an extended contact map.

Fig. S12. The effect of violating critical model assumptions.

Table S1. Rule-based model of the JNK network in neuroblastoma.

Table S2. Parameter estimates.

Table S3. Patient metadata.

Table S4. mRNA sequencing data.

Table S5. RPPA data.

References (6371)


Acknowledgments: We would like to acknowledge A. T. Look (Dana-Farber Cancer Institute, Harvard Medical School) for providing the transgenic Tg[dβh:MYCN-EGFP] zebrafish; W. Hahn and D. Root for providing the pDONR223-MKK7 construct; and W. Gallagher (Conway Institute, University College Dublin) for providing the pLKO control shRNA and pLKO ZAK shRNA. Data used in this publication were generated by the Clinical Proteomic Tumor Analysis Consortium (National Cancer Institute, NIH). Funding: The research leading to the results presented in this article was funded by the European Union Seventh Framework Programme (FP7/2007- 2013) ASSET project under grant agreement FP7- HEALTH-2010-259348, by Science Foundation Ireland under grant 06/CE/B1129, and by Breast Cancer Campaign UK under grant 2013MaySP024. N.R. was supported by the European Union Seventh Framework Programme (FP7/2007-2013) PRIMES project under grant agreement FP-HEALTH-2011-278568. M.F. was supported by German Cancer Aid (grant no. 110122), German Ministry of Science and Education (BMBF) as part of the e:Med initiative (grant no. 01ZX1303A and 01ZX1307D), Fördergesellschaft Kinderkrebs-Neuroblastom-Forschung e.V. D.R.C. was supported by Science Foundation Ireland under grant 11/SIRG/B2157 and by Cancer Institute NSW under grant 13/FRL/1-02. Author contributions: D.F. conceived the study, wrote and edited the manuscript, assisted with experimental design, and performed all mathematical modeling and survival analysis. M.H. performed zebrafish experiments. M.H. and J.F.H. performed Western blotting and coimmunoprecipitation experiments for Akt crosstalk and JNK feedback analysis, respectively. S.P.K., J.F.H., N.R., A.G.M., and R.P. performed all cloning, mutagenesis, and DNA purification. D.D. performed the mRNA purification and gene expression analysis. M.F. facilitated the collection of patient samples and data. F.W. supervised D.D. and facilitated the collection of patient samples and data. W.K. supervised M.H., N.R., A.G.M., and R.P.; wrote and edited the manuscript; and assisted with experimental design. B.N.K. supervised D.F., edited the manuscript, and assisted with experimental design. D.R.C. supervised S.P.K. and J.F.H, conceived the study; wrote and edited the manuscript; and performed experiments relating to JNK network activation, feedback analysis, and apoptosis. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All raw and normalized microarray data are available at the ArrayExpress database (; accession: E-TABM-38, E-MTAB-161, E-MTAB-179, and E-MTAB-1781). The Tg[dβh:MYCN-EGFP] zebrafish require a material transfer agreement from Dana-Farber Cancer Institute, Harvard Medical School.
View Abstract

Stay Connected to Science Signaling

Navigate This Article