Research ArticleCancer

Integrating network reconstruction with mechanistic modeling to predict cancer therapies

See allHide authors and affiliations

Science Signaling  22 Nov 2016:
Vol. 9, Issue 455, pp. ra114
DOI: 10.1126/scisignal.aae0535

Translating genomic mutations into drug treatments

A holy grail of systems biology is having sufficient information to use mathematical models to predict therapeutic strategies for patients. Halasz et al. developed a computational approach that used published information and newly acquired experimental data to construct signaling network models that could be tested to discover new potential therapeutic strategies for cancer. They applied this approach to colorectal cancer cell lines and identified a specific network feedback event present in a subset of the cells that was predicted to cause resistance to drugs that target the growth factor receptor EGFR. Testing this prediction in a zebrafish tumor migration model was consistent with the predictions. Thus, modeling patient cell signaling data may eventually aid in identifying the best personalized treatment.

Abstract

Signal transduction networks are often rewired in cancer cells. Identifying these alterations will enable more effective cancer treatment. We developed a computational framework that can identify, reconstruct, and mechanistically model these rewired networks from noisy and incomplete perturbation response data and then predict potential targets for intervention. As a proof of principle, we analyzed a perturbation data set targeting epidermal growth factor receptor (EGFR) and insulin-like growth factor 1 receptor (IGF1R) pathways in a panel of colorectal cancer cells. Our computational approach predicted cell line–specific network rewiring. In particular, feedback inhibition of insulin receptor substrate 1 (IRS1) by the kinase p70S6K was predicted to confer resistance to EGFR inhibition, suggesting that disrupting this feedback may restore sensitivity to EGFR inhibitors in colorectal cancer cells. We experimentally validated this prediction with colorectal cancer cell lines in culture and in a zebrafish (Danio rerio) xenograft model.

INTRODUCTION

Signal transduction networks carry extracellular signals from the cell membrane to the nucleus and internal signals from organelles or related to cell stress or metabolic status through finely controlled networks of protein interactions (1). Many diseases, such as cancer, are consequences of harmful changes in these networks, caused by oncogenic mutations or aberrant expression of genes encoding the participant proteins (1). Understanding the functional consequences of these alterations is key to effective treatment (1). Systems biologists have developed an array of computational tools to gain understanding into the network rewiring and functional consequences of such rewiring. Although mechanistic models of network are useful tools (2–9), their construction requires comprehensive prior information about the components and wiring structure (that is, the topology) of the network being modeled (10, 11) and plenty of experimental data for model calibration. Their application to the signaling networks in cancer cells is challenging, because these networks are often rewired with often unknown effects on the network (1).

An alternative strategy is to use experimental data and data-driven modeling instead of prior knowledge (12–24). However, such reconstructed networks, although useful in visualizing static network topologies, have limited use in predicting functional properties, such as the temporal and spatial dynamics of signaling networks, the consequences of sequential or pulsed drug treatments (22), the phenotypic differences after treatments with different drugs that target the same molecules using different mechanisms (22), and the acquired resistance after treatment. Most of these functional properties and perturbation outcomes can be predicted using mechanistic models. However, building mechanistic models based on reconstructed networks is not trivial, because the existing network reconstruction methods have limitations. For instance, some algorithms cannot detect feedback loops (15, 16), some cannot determine the direction of protein interactions (23), and some can only indicate the presence and absence of interactions (15–18) but cannot quantify the influence of one protein on another (19–22). Noise in measurements and suboptimal experimental designs affect the reliability and accuracy of the reconstructed networks (25).

Here, we present a two-stage computational platform that enables the construction of mechanistic models based on reconstructed models of signal transduction networks. The first stage consists of a network reconstruction algorithm that can reconstruct cell- and tissue-specific models of networks by combining prior knowledge of generic signal transduction network topologies with cell- and tissue-specific experimental data. The reconstructed models reveal the cell- and tissue-specific topology of the network and quantify the strengths of its interactions. In the second stage, an ordinary differential equation (ODE)–based mechanistic model, encompassing all previously known and newly inferred interactions of the network, is developed. The parameters of this model are then calibrated to ensure that the interaction strengths between different components of the mechanistic model closely resemble those inferred by the network reconstruction algorithm. The calibrated ODE model is then used to predict the effects of therapeutic intervention that targets the networks in different cells or tissues.

We tested our method on a publicly available data set (22) that contains quantitative measurements of different components of the epidermal growth factor receptor (EGFR) and insulin-like growth factor 1 receptor (IGF1R) signaling networks in a panel of colorectal cancer (CRC) cells. The first stage of network reconstruction uncovered differences in the topologies of the EGF and IGF1 signaling pathways of different CRC cell lines; some of these network differences were associated with EGFR inhibitor resistance in subsets of CRC cells. In the second stage, we constructed a mechanistic model of the EGF and IGF1 signaling networks for the EGFR inhibitor–resistant CRC cell line HCT116. The model correctly predicted the responses of the EGFR and IGF1R networks to different individual and combinatorial perturbations, as well as predicted a mechanism to overcome resistance to EGFR inhibition. We validated these predictions in cultured CRC cells and in zebrafish xenograft models.

RESULTS

Computationally integrating network reconstruction methods with mechanistic modeling of signaling networks

Our computational framework (Fig. 1A) consists of two stages: (i) network reconstruction (Fig. 1B) using a Bayesian reformulation of modular response analysis (BMRA) (19, 20) and (ii) mechanistic model development and calibration using a Bayesian mechanistic modeling (BMM) algorithm (Fig. 1C). BMRA combines prior literature knowledge with quantitative noisy experimental data to generate and statistically assess potential network models (topology and interaction strengths) of signaling networks (see note S1 and Materials and Methods for details). We compared the accuracy and robustness of BMRA with several existing network reconstruction algorithms. We benchmarked its performance using five data sets published by the DREAM (challenge 4) consortium (http://dreamchallenges.org/) (see note S2 for details). When applied without the inclusion of any “prior” knowledge of the network topology, BMRA ranked second, third, and fourth (of 29) on each of the five data sets, and with the inclusion of noisy prior knowledge, BMRA performed better than all of the original 29 algorithms in the challenge (Table 1).

Fig. 1 Schematic outline of the computational analysis pipeline and its components.

(A) Workflow of the integrated analysis platform. (B) Workflow of the BMRA algorithm, which reconstructs network models of signal transduction networks from perturbation data sets. (C) Workflow of the BMM algorithm, which reconstructs ODE models of signal transduction networks using reconstructed network models.

Table 1

Performances of BMRA and the top five performers of the DREAM challenge.

View this table:

We also systematically evaluated the sensitivity of BMRA to noise and the lack of experimental perturbation data using simulated data of the mitogen-activated protein kinase (MAPK) pathway based on a previously published mathematical model (18) of the pathway. Using the model, we simulated 2300 steady-state perturbation response data sets, 100 replicate data sets in 23 categories, containing different levels of noise and data incompleteness (see note S3 for details). We compared the performance of BMRA on these data sets with other methods based on Lasso regression (26–29), linear matrix inequality (30–32), and MRA by Klinger et al. (22), the latter of which was originally used to analyze the same experimental data set as we use here. BMRA outperformed the other algorithms in most cases except at high levels of noise, where all methods collapsed (Fig. 2, A to C). These results indicate BMRA’s higher accuracy and robustness against noise, low numbers of replicates, or sparseness of perturbations in experimental data.

Fig. 2 Sensitivity of BMRA to noise and data incompleteness in comparison to other algorithms.

(A) Performances of BMRA, LMML, pLasso, and CORENet in terms of area under the receiver operating characteristic (AUROC) curve and area under the precision recall (AUPR) curve values at different levels of noise. (B) Performances of BMRA, LMML, pLasso, and CORENet in terms of AUROC and AUPR values for decreasing numbers of replicates. (C) Performances of BMRA, LMML, pLasso, and CORENet in terms of AUROC and AUPR values for decreasing numbers of perturbation experiments.

The next step is to “feed” the BMRA-reconstructed network topology model into the BMM algorithm to develop predictive mechanistic network models (Fig. 1C and note S4). Any interactions that were inferred from the BMRA and not yet supported in the published literature should be validated experimentally, and then an ODE-based mathematical model can be developed for the final validated topological model. The mathematical model formulates the rate of change in the concentration of each pathway component by a parametric function dependent on its own concentration and that of its regulators. The parameters of these rate functions are then calibrated using a Bayesian parameter calibration algorithm.

Reconstructing EGFR and IGF1R networks in a panel of CRC cells

We used the computational pipeline (see note S5 for details) to analyze a publicly available data set containing perturbation responses of EGFR and IGF1R networks in six genetically diverse CRC cell lines: HCT116, HT29, LIM1215, SW403, SW480, and RKO (22). This data set comprises measurements of phosphorylation changes in eight key signaling proteins [AKT Ser473, extracellular signal–regulated kinase 2 (ERK2) Thr185/Tyr187, MAPK kinase 1 (MEK1) Ser217/Ser221, p70S6K Thr421/Ser424, IGF1R Tyr1131, glycogen synthase kinase 3 α/β (GSK3α/β) Ser21/Ser9, IκBα Ser32/Ser36, and insulin receptor substrate 1 (IRS1) Ser636/Ser639] in response to individual treatment by four inhibitor drugs [targeting the kinases MEK, phosphatidylinositol 3-kinase (PI3K), IκB kinase (IKK), or GSK3α/β] followed by stimulation with transforming growth factor–α (TGFα; a ligand of EGFR) and IGF1. Because the RKO cells showed very little response to the drugs, we excluded them from our analysis.

Despite the size of the data set, it is still incomplete and insufficient for most network reconstruction algorithms (12–17, 19–22). Because of the technology used to acquire the data, the measurements have a lot of noise (33). However, the data set probes the EGFR and IGFR pathways, which are well described in the existing literature. Thus, this data set provides a suitable test bed for our computational pipeline, which is designed to integrate prior knowledge with noisy and insufficient data to reconstruct a predictive model of the network. The data set enables the investigation of a pressing issue in cancer biology, that is, the problem of EGFR inhibitor resistance in CRC. EGFR inhibition is one of the mainstays in the therapy of metastatic CRC (34). Unfortunately, about 30 to 50% of CRC patients have a mutated KRAS gene that renders them resistant to EGFR-targeted therapies (35). In addition, about half of the patients with wild-type KRAS will develop resistance to EGFR inhibition (36). Several studies indicate that the extensive cross-talk between EGFR and IGF1R might contribute to the acquired resistance to EGFR inhibition (37). Therefore, developing predictive models of these pathways in EGFR inhibitor–resistant CRC cells (such as HCT116 and SW480) may provide valuable insight into the resistance mechanisms of CRC.

First, we developed a literature-based generic signaling network for EGFR and IGF1R pathways consisting of all known interactions (note S6) between the measured proteins across a diverse range of cell types. The generic network model represents our prior knowledge of EGFR and IGF1R pathways and was used to formulate the prior distribution for the BMRA algorithm. With the prior distributions, we used BMRA to reconstruct separate models for networks activated by TGFα or IGF1 in five CRC cell lines from the corresponding perturbation data sets. This resulted in 10 network models, corresponding to the five CRC cell lines each stimulated by TGFα (Fig. 3A) or IGF1 (Fig. 3B).

Fig. 3 Network reconstruction of the EGFR and IGF1R pathways in colon cancer cells.

(A and B) Reconstructed network models of cells stimulated with TGFα or IGF1. The left graphs show the predicted interaction strengths, with red bars indicating positive interaction strengths and blue bars indicating negative interaction strengths. The x axis represents interactions that occur in the TGFα- or IGF1-stimulated cells, the y axis represents cell lines, and the z axis represents interaction strengths. The height of each bar and the associated error bar represent the absolute mean and SD of the corresponding interaction strength. The width of each bar represents the corresponding interaction probability. The right network diagrams show the differences in the topologies of the TGFα- or IGF1-stimulated pathways in different CRC lines. EGFR is not part of the reconstructed networks, because it was not measured in the perturbation data set (22). We included it in these diagrams for convenience of visualization. The squares and circles associated with each connection represent the positive (+ve) or negative (−ve) interaction strengths. The size of the circles or squares indicates the interaction strengths; the cell lines are color-coded.

In the reconstructed networks, each interaction is annotated with cell-specific symbol for the probability distribution of its strength, which includes whether the interaction has a positive or negative effect, and the probability of its existence (Fig. 3, A and B). The reconstructed models predicted several differences in the signal integration, interpathway links, and feedback mechanisms of the EGFR and IGF1R pathways across the different CRC cell lines. For instance, the networks indicated that ERK inhibits AKT phosphorylation in LIM1215 and SW403 cells stimulated with TGFα, but this inhibition was predicted to be much weaker in HCT116 and SW480 cells. Furthermore, the networks predicted that ERK stimulates AKT phosphorylation in HT29 cells in response to TGFα. The networks of IGF1-stimulated cells indicated that ERK inhibits AKT phosphorylation in all cells, with the weakest effect predicted for SW480 cells. In IGF1- or TGFα-stimulated cells, the networks predicted that AKT promotes MEK phosphorylation in all cells except the SW480 cells, in which the relationship from AKT to MEK is negative in IGF1-stimulated cells. The networks also predicted differences in signal integration resulting from the actions of ERK and AKT on GSK3β and p70S6K phosphorylation. The HCT116 cell networks predicted that GSK3β is regulated by both ERK and AKT, but the LIM1215 cell networks indicated mostly ERK as promoting GSK3β phosphorylation, and the HT29 cell networks predicted an inhibitory effect of AKT on GSK3β phosphorylation in response to TGFα and a relatively small positive effect in response to IGF1. Similarly, in SW403 cells, the networks predicted differential AKT-mediated effects on GSK3β in response to IGF1 or TGFα: IGF1 promoted a positive and predominant regulation of GSK3β phosphorylation by AKT, and TGFα produced an inhibitory effect of AKT on GSK3β phosphorylation. The SW403 cells were again predicted to be different from the other cell lines with regard to regulation of p70S6K phosphorylation, which is stimulated by ERK in all cell lines and also by AKT in all except the SW403 cells.

Feedback loops were also predicted to operate differentially between cell lines. For instance, p70S6K can function as part of a negative feedback loop that inhibits IRS1 by phosphorylating its inhibitory sites (Ser636 and Ser639). This negative feedback loop was strongest in the models of SW480 cells and was not present in the models of LIM1215 cells exposed to either TGFα or IGF1. Both HCT116 and SW480 are resistant to EGFR inhibitor drugs (38), and negative feedback loops may play important roles in drug resistance (39). Therefore, we chose this feedback loop for further exploration and performed biochemical experiments to (i) validate its existence and (ii) determine whether it is associated with any of the commonly occurring oncogenic mutations.

Differential feedback inhibition of IRS1 by p70S6K in CRC lines

We chose four cell lines for the biochemical experiments: HCT116 and SW480 cells, which have oncogenic KRAS mutations; HKE3 cells, which are genetically identical to HCT116 with the exception that the mutant KRAS allele was knocked out, leaving the cells with a single wild-type KRAS allele (40); and HT29 cells, which have a BRAFV600E mutation. We investigated the strength of the p70S6K to IRS1 feedback loop by inhibiting p70S6K phosphorylation or reducing p70S6K abundance and subsequently measuring IRS1S636/S639 phosphorylation 0, 10, and 60 min after TGFα stimulation. We inhibited p70S6K phosphorylation using a small-molecule inhibitor of AKT, which directly phosphorylates p70S6K, and small interfering RNA (siRNA) that reduces p70S6K abundance. The AKT inhibitor reduced TGFα-stimulated p70S6K phosphorylation in all four cell lines (Fig. 4A and fig. S2), although the effect was less in the SW480 cells. Only in the HCT116 and SW480 cell lines did AKT inhibition significantly reduce IRS1S636/S639 phosphorylation. The siRNA-mediated p70S6K knockdown was variably effective across the cell lines, with the HCT116 and HKE3 cells showing the most effective reduction in p70S6K abundance. Despite the differences in knockdown efficiency, the p70S6K-targeted siRNA caused a significant reduction in IRS1S636/S639 phosphorylation in HCT116 and SW480 cells, but not in HKE3 and HT29 cells, 60 min after TGFα stimulation (Fig. 4B and fig. S3). Thus, the observed differences in IRS1S636/S639 phosphorylation did not correlate with efficiency of knockdown and were more likely to represent cell line differences in the strength of the feedback loop.

Fig. 4 p70S6K-mediated feedback inhibition of IRS1 in SW480, HCT116, HKE3, and HT29 cells.

(A) Effects of AKT inhibition on p70S6KT421/S424 and IRS1S636/S639 phosphorylation in SW480, HCT116, HKE3, and HT29 cells. Starved cells were treated with AKT inhibitor VIII (10 μM) for 1 hour and then stimulated with TGFα (100 nM) for 0, 10, and 60 min. Phospho-AKTS473, phospho-IRS1S636/S639, total IRS1, phospho-p70S6KT421/S424, and total p70S6K were measured by Western blotting. (B) Effects of p70S6K knockdown (KD) on IRS1 phosphorylation in SW480, HCT116, HKE3, and HT29 cells. Cells were transfected with siRNA against p70S6K. Twenty-four hours later, cells were serum-starved for 4 hours, treated with TGFα, and analyzed as above. Blots were quantified using ImageJ. In (B) and (C), the phosphorylated proteins were normalized to the respective total proteins, and the normalized levels were scaled between 0 and 1 and plotted. Error bars were calculated from n= 3 independent experiments (figs. S2 and S3). P values were calculated using Kruskal-Wallis test. (C) Graphical summary of the results. The squares and circles represent the positive (+Ve) or negative (−Ve) interaction strengths, respectively. GAPDH, glyceraldehyde-3-phosphate dehydrogenase.

These observations suggested that the p70S6K to IRS1 negative feedback is prominent in HCT116 and SW480 cells, but it is weaker or not present in HT29 and HKE3 (Fig. 4C). The differences in this feedback between HCT116, SW480, and HT29 cells agree with the results of the networks resulting from the BMRA algorithm (Fig. 3A). Additionally, the differences in this feedback between HCT116 and HKE3 cells, which differ only by the presence or absence of the mutated KRAS allele, suggested that the KRAS mutation contributes to the strength of this negative feedback loop. However, the BMRA algorithm did not predict that this feedback would be strong in SW403 cells, which also have the KRAS mutation. This may be due to noise in the perturbation data set that prevented the BMRA algorithm from detecting this feedback in SW403 cells. Alternatively, the KRAS mutation alone may not suffice to enable this feedback loop, which may depend on other genotypic features that are common to HCT116 and SW480 cells but are absent in SW403 cells. However, the data indicated that, in the context of these CRCs, the KRAS mutation is necessary but may not be sufficient for the p70S6K to IRS1 negative feedback loop.

Both HCT116 and SW480 contain the negative feedback loop and are resistant to EGFR inhibitor drugs (38). However, the network models reconstructed by the BMRA algorithm cannot provide a mechanistic explanation of whether and how this feedback loop, in association with the KRAS mutation, causes resistance to EGFR inhibitors. The analysis so far also does not enable us to predict whether targeting this feedback will sensitize the cells to subsequent treatment with EGFR inhibitors. Therefore, we used the BMM algorithm to derive a mechanistic model of the EGFR and IGF1R pathways in HCT116 cells and used this mechanistic model to address these questions.

Mechanistic modeling of EGFR and IGF1R networks in HCT116 cells using the BMM algorithm

The main purpose that we had for developing the mechanistic model was to study the effects of ligands (TGFα for the EGFR pathway and IGF1 for the IGF1R pathway), EGFR inhibitors, and perturbations to p70S6K to IRS1 negative feedback loop on the amount of active ERK and AKT in HCT116 cells. We developed an ODE model of the generic EGFR and IGF1R pathways using the generic network (Fig. 5A). We could use the generic network, because BMRA did not infer any previously unknown interactions (Fig. 3).

Fig. 5 Mechanistic modeling of EGFR and IGF1R pathways in HCT116 cells.

(A) Extending the generic signaling network for mathematical modeling. The indexed gray dashed lines in the left panel represent the interactions that were expanded for ODE modeling. The blue lines in the right panel represent the corresponding (indicated by the index) expanded interactions. The red lines in the right panel represent newly added interactions. The light gray nodes and lines in the right panel represent the interactions that were excluded from the ODE modeling. (B) Mathematical model of the generic EGFR and IGF1R pathways. Each protein has an inactive and active state denoted by the suffix “a.” Proteins that are phosphorylated at inhibitory sites also have inhibited states denoted by the suffix “i.” MM(kx, Kx) represents the Michaelis-Menten kinetics with rate and dissociation constants kx and Kx, respectively; MM1(KMx, VMx) represents the Henri-Michaelis-Menten kinetics with the Michaelis constant KMx and the maximum velocity constant VMx; HM(kx, Kx, Kx1, β) represents the general hyperbolic modifier kinetics (29) with the rate constant kx, dissociation constant Kx, activation or inactivation constant Kx1, and modification constant β. The molecules in red have activating mutations in HCT116 cells. The model was found to be sensitive to the parameters highlighted in yellow. (C) Prior and posterior distributions of the sensitive parameters. The x axis of each subplot represents the values of the corresponding parameter, and the y axis represents the posterior probabilities associated with these values. The distributions of prior values are represented by gray curves.

The generic model (Fig. 5A) only contains a single receptor (IGFR) and lacks RAS, which we determined was important for the negative feedback loop and is present in an oncogenic form in HCT116 cells, and the kinase RAF, which is activated by RAS and signals to MEK and AKT (3, 22, 41). The generic model also lacks PI3K, which is mutated in HCT116 cells (3, 41). Therefore, we had to modify the generic network to include both receptor inputs EGFR (the target of EGFR inhibitors) and IGFR, RAS (representing KRAS and connecting EGFR to RAF), RAF (connecting RAS and AKT to MEK), and PI3K (connecting IRS to downstream AKT). To reduce the complexity of the model, we excluded GSK3β and IκB, which are downstream of ERK and AKT and have limited influence on the activities of either ERK or AKT. The resulting model (Fig. 5B) contains 10 proteins—EGFR, IGFR, RAS, RAF, MEK, ERK, IRS, PI3K, AKT, and p70S6K—and the relationships connecting them.

We collected information about the mechanism of protein activation or inhibition from the literature (3–6, 8–10, 24, 34). Protein activation or inhibition events are often multistage processes that can involve multisite phosphorylation, conformation change, dimerization, or interaction with other partners, among others (3–6, 8–10, 24, 34). To reduce model complexity, we simplified these processes to single partial steps (3) that transform a protein from its inactive to active or inhibited form or vice versa. We modeled most protein activation and inhibition processes using Michaelis-Menten kinetics, except for the modulation of AKT phosphorylation by ERK, which is indirect and occurs through unknown components. We modeled the ERK-mediated regulation of AKT as a hyperbolic modifier function (42).

The resulting ODE model had 66 parameters that were initially unknown (Fig. 5B). The BMM algorithm identified 17 parameters (highlighted yellow in Fig. 5B) that were predicted to have the most impact on the model output (note S6). We used the BMM algorithm to estimate HCT116-specific posterior distributions of these parameters while keeping the remaining parameters at fixed values (note S6). BMM iteratively compared the simulated interaction strengths among IGFR, IRS, AKT, MEK, ERK, and p70S6K with those estimated by the BMRA algorithm (Fig. 3, A and B) for parameter estimation. We excluded the interaction strengths among GSK3β, IκB, AKT, and ERK from the BMM algorithm, because these proteins are not in the ODE model. The estimated posterior distributions (Fig. 5C) reflect well-known characteristics of EGFR and IGF1R pathways in HCT116 cells. For instance, the activating KRAS mutation in HCT116 cells was well described by the posterior distributions of the RAS activation and inactivation rates, which indicated that RAS is activated at a much higher rate (kras1 = 1.35 ± 0.36 min− 1 and kras2 = 3.02 ± 1.5 min−1; see Fig. 5C) than when it is inactivated [VMras = 0.012 ± 0.18 arbitrary unit (AU) min−1]. Additionally, the inferred posterior distributions of many model parameters have more than one mode or peak and, in most cases, have two dominant peaks (bimodal) (Fig. 5C), which is consistent with what is known about how the activation and inactivation rates of the components vary depending on the input signal and the stimulation of the EGFR or IGFR (3, 9, 41).

Initial validation of the mechanistic model

We performed model simulation using an ensemble of ODE models, each of which was fitted with parameter values sampled from the corresponding posterior distribution (note S6). We then used the average and SD of the ensemble simulations to represent simulation results and the corresponding confidence intervals. First, we verified that the interaction strengths simulated by the calibrated model resembled those estimated by the BMRA algorithm (fig. S1, A to C). Subsequently, we simulated the dynamics of the amount of active ERK and AKT in response to different amounts of ligands. The model simulations suggested that ERK and AKT were transiently activated at low levels of ligand stimulation but then quickly achieved sustained activation at higher levels (Fig. 6A). To validate these predictions, we tested the dynamics of AKT and ERK phosphorylation in serum-deprived HCT116 cells in response to high and low concentrations of EGF and insulin. EGF and insulin activate the same receptors as TGFα and IGF1, but EGF and insulin are more commonly used in literature for studying ERK and AKT activities. By using EGF and insulin for these experiments, we could use existing information from literature to determine different aspects for the validation experiments, such as high and low ligand concentrations and optimal time points for measurement. The measured dynamics of ERK and AKT in response to the two different concentrations of EGF were similar to those predicted by the model simulations—transient and sustained activation at low and high ligand concentrations, respectively (Fig. 6B and fig. S4). The response to insulin was not as consistent with the simulations as the EGF response was (Fig. 6C and fig. S4). In agreement with the model simulation, the cells exhibited sustained phosphorylation of AKT at high insulin concentration. However, at low insulin concentration, simulations indicated an initial peak in AKT activity at ~10 min, which then dipped at ~40 min before mostly recovering at 60 min (Fig. 6A). This transient reduction in AKT phosphorylation was not captured in the experimental results, which were obtained at 0, 10, and 60 min (Fig. 6C and fig. S4). Also, the cells had high basal ERK phosphorylation even before insulin stimulation. This was an unexpected observation and did not agree with simulations, because in the simulation study, we assumed that starved cells start from a basal state lacking phosphorylation (2–9).

Fig. 6 Biochemical validation of the mechanistic ODE model.

(A) Simulated dynamics of ERK and AKT activation for different levels of EGF/TGFα and insulin/IGF1 stimulations, respectively. The x axis represents time in minutes, the y axis represents EGF or insulin levels in AUs, and the z axis represents the concentrations of activated ERK (A and C) and AKT (B and D) in AUs. Solid lines represent average activities, and the shaded areas correspond to SDs. (B and C) Phosphorylation levels of ERK1/2T202/Y204 and AKTS473 at 0, 10, and 60 min were determined by Western blotting after stimulation of serum-starved HCT116 cells with EGF (20 or 100 ng/ml) (B) or insulin (1 or 10 μg/ml) (C). Blots were quantified using ImageJ. The phosphorylation levels of ERK1/2T202/Y204 and AKTS473 were then normalized to the loading control GAPDH, scaled between 0 and 1, and plotted. The data represent n = 3 independent experiments (fig. S4), with error bars indicating SD. P values were calculated using Kruskal-Wallis test.

Despite some discrepancies between the simulation and experimental data, which can be attributed to unexpected basal phosphorylation and low temporal resolution in the experimental data, the model did match ERK and AKT dynamics in many respects in response to ligand stimulations. Therefore, we thought that the model would be useful to explore the role of the p70S6K to IRS1 negative feedback loop in the resistance of HCT116 cells to EGFR inhibitors.

Sensitizing the model to EGFR inhibitors by simulating p70S6K knockdown

Model simulations indicated that knocking down p70S6K would not affect ERK phosphorylation but would increase AKT phosphorylation (Fig. 7A) in HCT116 cells stimulated through the EGFR pathway. To verify this prediction, we knocked down p70S6K in the serum-starved HCT116 cells and exposed the cells to TGFα (Fig. 7B and fig. S5). Although both EGF and TGFα were previously shown to strongly stimulate p70S6K phosphorylation (43, 44), we used TGFα in these experiments because we detected the p70S6K- and IRS1-mediated feedback loop in TGFα-stimulated networks (Fig. 3). p70S6K knockdown increased the amount of TGFα-stimulated AKT phosphorylation without causing a significant change in ERK phosphorylation (Fig. 7B and fig. S5), suggesting that AKT activity is affected by the p70S6K to IRS1 negative feedback loop. Hence, we analyzed AKT activity in serum-grown HCT116 cells using inhibitors or knockdown to perturb the cells without exogenous ligands in addition to the ones in the serum.

Fig. 7 The role of the p70S6K-mediated negative feedback to IRS1 in EGFR inhibitor resistance of HCT116 cells.

(A) Simulated active ERK (blue) and AKT (red) concentrations after p70S6K knockdown and different levels of EGF or TGFα (EGF/TGFα) stimulations. Solid lines and shaded areas represent means and SDs. (B) Effects of p70S6K knockdown on AKTS473 and ERK1/2T202/Y204 phosphorylation. HCT116 cells were transfected with nontargeting siRNAs (CTRL) or siRNAs against p70S6K (KD) and grown in serum (see fig. S2 for knockdown efficiency) for 24 hours. Then, starved (4 hours) cells were stimulated with TGFα (100 nM) for the indicated time points. Protein levels were quantified and normalized as above. Data are representative of n = 3 independent experiments (fig. S6). P values were calculated using Kruskal-Wallis test. (C) Simulated effects of p70S6K knockdown and EGFR inhibition on AKT activation in serum-grown HCT116 cells. The starting amounts of ligand were set to 0.01 AU. (D and E) Effects of p70S6K knockdown and EGFR inhibition on AKTS473 phosphorylation in serum-grown HCT116 cells. Control and knockdown HCT116 cells were treated with the EGFR inhibitor BIBX1382 (BIBX; 5 μM). Phosphorylated AKTS473 levels were measured at the indicated time points after treatment and normalized by total AKT levels. Data are from n = 3 independent experiments (fig. S5). P values were calculated using Kruskal-Wallis test. A representative blot is shown in (E). (F) Simulation of AKT activity for different p70S6K knockdown efficiencies and EGFR inhibitor strengths. (G) Phosphorylated AKT S473 levels in control and knockdown HCT116 cells with and without BIBX1382 treatment. Data represent n = 3 independent experiments. P values were calculated using Kruskal-Wallis test. (H) Feedback loops that control AKT phosphorylation. Positive feedback loops (highlighted in blue) controlling AKT phosphorylation and negative feedback loops (highlighted in red) controlling AKT phosphorylation.

In the simulations of ligand-stimulated responses, the system starts from a resting point with no active AKT or ERK. When we simulate the effects of p70S6K knockdown or EGFR inhibition in the presence of serum, the system is not starting from baseline of no activity. First, we simulated active AKT concentration in response to either p70S6K knockdown or EGFR inhibition. As in the case of ligand-mediated stimulation of the EGFR pathway with either TGFα or EGF, simulating p70S6K knockdown produced an increase in the amount of active AKT; simulating EGFR inhibition with a reversible and selective inhibitor showed a small reduction in the amount of active AKT at ~12 min before recovering to its initial amount at ~27 min (Fig. 7C). Analysis of AKT phosphorylation in HCT116 cells exposed to the reversible and selective EGFR inhibitor BIBX1382 or in which p70S6K was knocked down showed profiles consistent with the simulation results (Fig. 7, D and E, and fig. S5). We then simulated the amount of active AKT in response to combinations of different p70S6K knockdown efficiencies and “doses” of EGFR inhibitor (note S6) in serum-grown HCT116 cells. The simulation suggested that a combination of partial p70S6K knockdown and EGFR inhibitor results in reduction in AKT phosphorylation compared to either individual treatments (Fig. 7F). This was confirmed in biochemical experiments (Fig. 7, D and G, and fig. S5). This effect on AKT activation is counterintuitive, because one might predict that loss of a negative feedback loop would enhance activity rather than suppress it.

The presence of multiple positive and negative feedback loops that affect AKT activity may explain this counterintuitive behavior of AKT in response to the combination of p70S6K knockdown and EGFR inhibition. AKT activity is tightly controlled by several negative and positive feedback loops (Fig. 7H), and some of these involve p70S6K. In the HCT116 cells, the negative feedback loop involving AKT and p70S6K is usually operational and causes an increase in AKT phosphorylation in response to p70S6K knockdown (Fig. 7, A to E). However, the positive feedback loops involve the RAF-MEK-ERK pathway, which becomes saturated in the presence of ligand, aided by mutant RAS, and is not affected by p70S6K knockdown (Fig. 7, A and B, and fig. S6). EGFR inhibition releases this pathway from this saturating effect, enabling the positive feedback loops involving p70S6K and AKT to have an effect. When both positive and negative feedbacks involving p70S6K and AKT are operational, AKT phosphorylation may either increase or decrease in response to p70S6K knockdown depending on which feedback loops have more impact on its phosphorylation. The positive feedback loops involving p70S6K are reinforced by several other positive feedback loops (Fig. 7H), involving the RAF-MEK-ERK pathway. Therefore, it is possible that, in the presence of EGFR inhibitors, the positive feedback loops have a stronger impact on AKT phosphorylation than the negative feedback loops, which leads to a decrease in AKT phosphorylation in response to p70S6K knockdown.

Synergistic effects of EGFR inhibition and p70S6K knockdown in culture and in vivo

Inhibition of AKT pathway increases apoptosis and decreases invasiveness of some cancer cells (45–50). Therefore, we tested the effects of combining EGFR inhibition with p70S6K knockdown on serum-grown HCT116 cells (Fig. 8A and fig. S7). The combination induced a significantly higher level of apoptosis in HCT116 cells compared to either individual treatment. The Bliss independence score of the combination, which is the log2 ratio of the expected to observed cell mortality rates, was found to be −0.217, suggesting that the combined treatment killed more cells than the individual treatments combined and indicating a synergistic effect of the combination on HCT116 cell death. HKE3 and HT29 cell lines, which our analysis indicated lack this feedback loop, exhibited increased apoptosis in response to EGFR inhibition, and the amount of toxicity was unchanged by p70S6K knockdown (Fig. 8A and figs. S8 and S9).

Fig. 8 Synergism between EGFR inhibition and p70S6K knockdown in cultured and xenografted HCT116 cells.

(A) Effects of p70S6K knockdown and EGFR inhibition on apoptosis of CRC cells. Left, HCT116; middle, HKE3; right, HT29. Control and knockdown cells were treated with 5 μM BIBX1382, and apoptosis was measured 48 hours later. Fold changes in apoptosis with respect to control cells are shown. Results represent four independent experiments (figs. S7 to S9). P values were calculated using Kruskal-Wallis test. (B) Workflow of the zebrafish xenograft–based model validation system. dpf, days post-fertilization. (C) Effect of lapatinib or control [dimethyl sulfoxide (DMSO)] treatment on zebrafish embryos co-injected with DiO-labeled control and DiI-labeled p70S6K knockdown HCT116 cells in the yolk sac. Fish were assessed for survival and dissemination of injected HCT116 cells (the number of fish with no cells detected outside the yolk sac = no dissemination or cells detected in the tail) 3 days after injection and treatment. (D) Quantitation of disseminated cells within the zebrafish after DMSO (n = 11 fish) or lapatinib (n = 9 fish) treatment (mean ± SEM). P values were calculated using Kruskal-Wallis test. (E) Microscopic images of disseminated DiO-labeled control (CTRL; green) and DiI-labeled p70S6K knockdown (KD; red) HCT116 cells in the zebrafish tail after DMSO or lapatinib treatment. The white bar represents 0.5 mm.

We also tested the effect of the EGFR inhibitor lapatinib on metastasis of control or p70S6K knockdown HCT116 cells in a zebrafish (Danio rerio) metastasis assay. We used lapatinib because we found that BIBX1382 was toxic to zebrafish embryos, which died shortly after exposure to BIBX1382. The dye-labeled HCT116 cells that migrated out of the yolk sac (the site of injection) to the tail represented metastasis. Lapatinib treatment of embryos injected with control HCT116 cells reduced the number of fish with cells in the tail (Fig. 8, B and C). Lapatinib treatment was even more effective when used on the embryos injected with p70S6K knockdown cells (Fig. 8D), consistent with the synergistic effect that we observed in the HCT116 cells in culture.

These results suggested that the negative feedback loop from p70S6K to IRS1 plays a role in defining the sensitivity of HCT116 cells to EGFR inhibition line and that inhibition of p70S6K could increase the cytotoxic effects of EGFR inhibition. This further strengthens the hypothesis that p70S6K-mediated negative feedback loop plays a role in the EGFR inhibitor resistance of HCT116 cells.

DISCUSSION

Here, we present a computational pipeline that reconstructs network models by combining prior knowledge with perturbation data sets and then develops predictive mechanistic models on the basis of the reconstructed networks. This method of developing predictive network models relaxes the dependence of such predictive models on prior knowledge and provides a data-driven approach. This approach also overcomes some limitations, for example, sensitivity to noise and data incompleteness, of other network reconstruction and modeling approaches. Furthermore, BMRA characterizes each interaction by the probability distribution of its strength and the probability of its existence. These attributes reveal critical details of each interaction: its strength (mean of its distribution), type [that is, whether the interaction activating or inhibiting (sign of the mean)], direction, confidence interval for its strength (SDs of the distribution), co-dependence between different interactions (covariance between interaction strengths), and the probability of whether it is genuine or an artefact of spurious noise. Such details are invaluable in building mechanistic models based on reconstructed networks. A benefit of our mechanistic modeling approach is that it does not require extensive multiconditional time course experiments for model calibration. It relies on the quantitative and qualitative properties of the reconstructed networks.

Nevertheless, our approach is not without limitations: It is applicable in cases where a reasonable amount of prior knowledge and experimental data are available, and the ODE model in our method is calibrated to exhibit similar interaction strengths as the probabilistic networks, which are estimated from steady-state perturbation data. Therefore, these models are more likely to correctly simulate steady-state behavior than pathway dynamics. This could be remedied by using additional time course data for calibration. Finally, our algorithm is computationally intensive and, in its current form, can be used to analyze networks containing a maximum of 25 nodes within reasonable time. For bigger networks, the number of parameters that need to be calibrated increases, which would require using specialized Markov chain Monte Carlo (MCMC) methods designed for exploring high-dimensional parameter spaces, such as Riemann manifold Langevin–adjusted Metropolis-Hastings methods (51). Despite these limitations, we demonstrated that the current computational approach is a useful tool for exploring the mechanistic roles of signaling networks in disease and drug resistance and for predicting new therapeutic targets.

We demonstrated the effectiveness of our method by revealing a mechanism of EGFR inhibitor resistance in HCT116 cells using a small set of highly incomplete and noisy data. Additionally, our computational pipeline pinpointed a potential therapeutic approach to overcome EGFR inhibitor resistance in colon cancer cells. We identified a role of the p70S6K to IRS1 negative feedback loop in EGFR inhibitor resistance of metastatic CRC lines. Here, we provided experimental support for the role of this feedback in EGFR inhibitor resistance and identified an association of this feedback loop with KRAS mutation in some colon cancer cells. This feedback loop has been suspected, but not supported experimentally, in EGFR inhibitor resistance of several breast, colon, lung, kidney, pancreas, and glioblastoma cancer cells (52, 53). Because inhibition/siRNA-mediated knockdown of p70S6K synergistically kills head neck and small cell lung cancer cells when combined with EGFR inhibitors (54, 55), this feedback loop may be an important target for many different cancers. Here, using multiple methods, we showed that p70S6K knockdown synergizes with subsequent EGFR inhibitor treatment to kill KRAS-mutated colon cancer cells, but not cells lacking this mutation. Additionally, this combination synergized in reducing metastasis of KRAS-mutated colon cancer cells in zebrafish xenografts. Therefore, our study contributed valuable insight for treating patients with KRAS-mutated, EGFR inhibitor–resistant, metastatic CRC.

MATERIALS AND METHODS

The BMRA algorithm

Consider an N node (ni, i = 1 … N) biochemical network. MRA (19) showed that following a perturbation (pk), the fractional change (Rik=d log(xiSS)dpk) in the steady-state concentration (xiSS) of a network node (ni) is linearly related to that (Rjk, j ≠ i) of the other nodes (nj, j ≠ i) through the corresponding interaction strengths (rij, j ≠ i).∑j=1,j≠inrijRjk=Rik;i≠k;i,k=1…N(1)

If node nj, j ≠ i directly influences node ni, then rij ≠ 0; otherwise, rij = 0. We modified this equation by introducing two more variables, a binary variable Aij in Eq. 1 to indicate whether nj directly influences ni (Aij = 1) or not (Aij = 0) and an error variable (ϵik) to account for the imbalance in Eq. 1 caused by noisy experimental data. The resulting modified MRA equation is shown below.∑j=1,j≠inAijrijRjk+ϵik=Rik;i≠k;k=1…N(2)

We used Bayesian statistics to infer the probability distributions of Aij, i ≠ j and rij, i ≠ j from experimentally measured R. First, we conceptually divided the N node network into N smaller subnetworks, each of which consists of a specific node (ni) and its regulators. The probability distributions of the network descriptors (Ai = {Aij, j ≠ i}, ri = {rij, j ≠ i}) of each of these subnetworks were inferred separately. We first formulated prior probability distributions of Ai and ri based on existing knowledge of the network topology. The binary vector Ai denotes the known regulators of node ni. The prior distribution (P(Ai)) of Ai was then formulated as follows:P(Ai)∝exp(−ψ×dH(Ai,Ai))(3)where ψ is a constant and dH(Ai,Ai) is the Hamming distance between Ai and Ai. P(Ai) favors the subnetworks that are in agreement with existing knowledge and penalizes those that are in disagreement. The extent of penalty depends on ψ with larger values, resulting in stiffer penalties. The prior distribution of ri is dependent on Ai and is denoted by P(ri|Ai). In the absence of a direct interaction from nj to ni (that is, Aij = 0), the corresponding interaction strength (rij) was assumed to have 0 value with probability 1, whereas the interaction strengths representing direct interactions (Aij = 1, j ≠ i) were assumed to have Gaussian priors (note S1) (18, 56). The above prior distributions (P(Ai), P(ri|Ai)) were then used to calculate the likelihood (P(R|ri, Ai)) of the global changes (R) given the network descriptors (Ai, ri) (note S1). Subsequently, we applied Bayes’ rule (57) to calculate the joint posterior distribution (P (Ai,ri |R)) for the network descriptors (ri , Ai). Using the chain rule, the above joint posterior P(Ai, ri |R)) is the product of the marginal posterior (P(Ai|R)) of the binary indicators (Ai) and the conditional posterior (P(ri |Ai, R)) of the interactions strengths (ri), that is, P(Ai, ri |R) = P(ri| Ai, R) × P(Ai|R). The marginal posterior (P(Ai|R)) is proportional to the product of a multivariate Student’s t (MVSt) and the prior distribution, whereas the conditional posterior (P(ri| Ai, R)) is proportional to the MVSt distribution (note S1). The constant of proportionalities of the above posterior could not be calculated analytically; therefore, we used an MCMC algorithm to approximate the true joint posterior P(Ai, ri |R). The MCMC algorithm starts with a random network, and each iteration (t) randomly proposes a new subnetwork (Aitn) of Ai, which differs from the current subnetwork by only one interaction (note S1) (15, 58). The proposed subnetwork was then accepted (Ait=Aitn) with probability α= min(1,P(Aitn|R)P(Ait|R)) (note S1) (59). Subsequently, we drew a sample (rit) from the conditional distributions P(ri|Ait,R) using a random vector generator that samples from custom MVSt distributions. The above two steps were repeated for a large number (105) of iterations, and after an initial burn-in period (5 × 104 iterations), the accepted network topologies and interaction strengths were recorded. The mean of an element Aij of Ai over the post–burn-in samples was then used as the posterior probability of the interaction from node nj to ni, that is, P(Aij|R)=15×104∑t=5×104+1105Aijt. On the other hand, the sample means, variance, and covariance of different interaction strengths (rij) across the post–burn-in samples were used as the posterior mean, variance, and covariance, respectively, of the strengths of the interaction from node nj to ni. A MATLAB implementation of the above (BMRA) algorithm can be found in the Supplementary Materials (code S1).

Calibrating the parameters of the ODE model

To calibrate the parameters (Θ) of the ODE model, we first identified those parameters (Θs ⊆ Θ) that had the maximum impact on the model output(s). Here, the model output (Mo) is defined as the area under the concentration versus time curve (xo(t)) of the output node (no) of the model, that is, Mo=∫t=0∞xo(t)dt. When a model has many output nodes (no1, no2 …), the sum of the areas of their concentration curves (xo1(t), xo2(t), …) is used as the model output (Mo=∫t=0∞(xo1(t)+xo2(t)+…)dt). If a small change (Δθk) in a parameter (θk) causes the model output to change by ΔMo, then the sensitivity (Sk) of the model to this parameter is given by Sk= log(|ΔMoΔθk|). These sensitivities quantify the impact of a model’s parameters on its output. Model sensitivities (Sk) were estimated for each model parameter (θk) for a large number (105) of randomly generated (note S6) parameter settings using Monte Carlo sampling (fig. SN6.1). The parameters that had low (≈ 0) average sensitivities were assigned fixed values (note S6), and those with higher average sensitivities (≫ 0) were calibrated using a Bayesian parameter estimation algorithm. First, we assigned a prior probability distribution (P(Θs))(note S6) to the sensitive parameters (Θs) and formulated a likelihood function (L(Θs)) for the same: L(Θs)∝∏l|Σrl|−12exp(−(rMl−μrl)TΣrl−1(rMl−μrl)). Here, rl represents the simulated interaction strengths in condition l (for example, ligand stimulation); μrl and Σrl are their means and covariance matrices as estimated by the BMRA algorithm from experimental data. We then used Bayes’ rule (57) to calculate the posterior distribution [P(Θs| μrl, Σrl, l = 1, 2 …)] of the sensitive parameters (Θs): P(Θs|μrl,Σrl,l=1,2…)∝L(Θs)P(Θs)˙ However, the posterior distribution could only be calculated up to a constant of proportionality. Therefore, MCMC sampling was used to estimate the true posteriors of these parameters. The MCMC algorithm starts with a random set of values for the sensitive parameters, and each iteration (t) proposes a new set of values (Θstn) based on their current values (Θst) using a Gaussian proposal distribution (q(Θstn|Θst)), which is defined in log space, that is, log(Θstn)∼N(log(Θst),VΘs), where VΘs is the covariance matrix of the proposal distribution (note S6). The new values (Θstn) are then accepted (Θst=Θstn) with probability α= min(1P(Θstn|μrl,Σrl,l=1,2…˙)P(Θst|μrl,Σrl,l=1,2…˙)) (note S6) (59). This process was repeated many times (2.5 × 105), and the samples drawn after a burn-in (105) period were used to estimate the posterior distribution of Θs. The variances (VΘs) of the proposal distribution were adaptively (note S6) tuned to improve mixing of the MCMC chains (60). For a robust exploration of the parameter space, we deployed multiple parallel Markov chains that operate at different “temperatures” (parallel tempering) (61). The chains that operate at high temperatures (temperature ≫ 1) take large steps and can explore a vast region in the parameter space relatively quickly. The chains that operate at low temperature (≈ 1) take small steps and perform local searches. These chains pass information among each other in a probabilistic manner, helping the Markov chains operating at low temperature to avoid getting stuck at local minima. The samples drawn by the chain operating at normal temperature (temperature = 1) were used to estimate the posterior distribution of Θs. A MATLAB implementation of the above calibration algorithm can be found in the Supplementary Materials (code S1).

Cell lines

The human colorectal adenocarcinoma cell line HT29 was acquired from the European Collection of Cell Cultures. The human colorectal carcinoma cell lines HCT116 and HKE3 were provided by S. Shirasawa (43). Cells were maintained in complete medium [high-glucose Dulbecco’s modified Eagle’s medium (DMEM) (Invitrogen Life Technologies), supplemented with 10% fetal bovine serum (FBS) (Invitrogen Life Technologies), penicillin (100 IU/ml), and streptomycin (100 μg/ml) (Invitrogen Life Technologies)] and incubated at 37°C in a humidified atmosphere of 5% CO2.

Reagents

The following inhibitors were used: PI3K inhibitor LY294002 (10 μM; Calbiochem), AKT inhibitor IV (10 μM; EMD Millipore), AKT inhibitor VIII (10 μM; Sigma-Aldrich), MEK inhibitor U0126 (10 μM; Promega), EGFR inhibitor BIBX1382 (10 μM; Calbiochem), and lapatinib (2.5 μM; Selleck Chemicals). Cells were treated with TGFα (100 ng/ml; PeproTech), EGF (100 ng/ml; PeproTech), or insulin (10 μg/ml; Sigma-Aldrich).

RNA interference

Validated siRNA to interfere with RPS6KB1 mRNA (sense, 5′-GGUUUUUCAAGUACGAAAAtt-3′; antisense, 5′-UUUUCGUACUUGAAAAACCtt-3′) and nontargeting control siRNA were predesigned and obtained from Ambion. Transfections were performed with Oligofectamine (Invitrogen Life Technologies) according to the manufacturer’s instructions. Briefly, diluted oligonucleotides were preincubated with Oligofectamine for 20 min at room temperature, and then the mixtures were added dropwise to 1.5 × 105 cells in Opti-MEM (Invitrogen Life Technologies). After 4 hours at 37°C, DMEM with 30% FBS was added to the cultures to achieve a final FBS of 10%. Cells were lysed after 24, 48, or 72 hours for Western blot analyses to test the efficiency of knockdown.

Quantitative immunoblotting

Cells were serum-starved for 4 hours before stimulation. Cells were lysed in 1× Lysis Buffer (Cell Signaling Technology Inc.) according to the manufacturer’s instructions. NuPAGE Bis-Tris gels (Novex, Life Technologies) were used to resolve proteins. Gels were transferred to polyvinylidene difluoride membrane at 30 V for 1 hour. Membranes were blocked with tris-buffered saline (TBS)–Tween (pH 7.4) containing 5% skim milk for 1 hour and incubated with 1:1000 diluted rabbit anti-human phospho-AKT (Ser473), AKT, phospho-p44/42 MAPK (Thr202/Tyr204), p44/42 MAPK, p70 S6 kinase, phospho-p70 S6 kinase (Thr389), IRS1, and phospho-IRS1 (Ser636/Ser639) antibodies or as control with mouse anti-human GAPDH (all from Cell Signaling Technology Inc.) in TBS-Tween with 5% bovine serum albumin overnight at 4°C. After washing, bound antibodies were detected with 1:5000 diluted horseradish peroxidase–conjugated anti-rabbit or anti-mouse immunoglobulin G (Cell Signaling Technology Inc.), followed by development with WesternBright chemiluminescent substrate (Advansta). Membranes were scanned using a high-sensitivity charge-coupled device camera (Chemi Image, Advanced Molecular Vision). Quantification of the bands was performed by using the ImageJ software.

Apoptosis assays

The FITC Annexin V Apoptosis Detection Kit (BD Biosciences) was used to detect apoptosis by flow cytometry. Cells were exposed to BIBX1382 treatment for 48 hours, harvested (including detached cells), and processed according to the manufacturer’s instructions. Briefly, cells were washed twice with phosphate-buffered saline (PBS) and resuspended in 1× annexin V binding buffer. Then, 1 × 105 cells were incubated with 5 μl of annexin V–fluorescein isothiocyanate and 5 μl of propidium iodide for 15 min at room temperature. Finally, 400 μl of 1× binding buffer was added. Labeled cells were analyzed with an Accuri C6 flow cytometer (Becton Dickinson Immunocytometry Systems) equipped with the Accuri C6 software program (Becton Dickinson) for data acquisition and analysis.

Alternatively, after the drug treatments, cells were fixed in ethanol and stained with propidium iodide. The sub-G1 population was determined using an Accuri C6 flow cytometer.

Zebrafish xenografts

Zebrafish (D. rerio) maintenance and experiments were approved by the University College Dublin Animal Research Ethics Committee. Wild-type (AB line) embryos were produced by natural spawning and kept on a 14-hour light/10-hour dark cycle at 28.5°C. Control or p70S6K knockdown HCT116 cells were labeled with DiO or DiI (Vybrant, Invitrogen Life Technologies) lipophilic fluorescent dyes, respectively, as described by Halasz et al. (62). Equal amounts of control and p70S6K knockdown HCT116 cells were mixed, resuspended in PBS, and kept on ice before microinjections. Two days after fertilization, dechorionated zebrafish embryos were anesthetized with 0.016% MS-222 (Sigma-Aldrich), and about 100 cells were transplanted into the yolk sac of embryos using a microinjector [PV830 Pneumatic PicoPump (World Precision Instruments) with M-152 manipulator (Narishige)] equipped with a glass capillary (World Precision Instruments). Then, embryos were incubated for 1 hour at 31°C and checked for cell presence. Fish with fluorescent cells outside the implantation area were excluded from further analysis. All other fish were incubated with DMSO or 2.5 μM lapatinib (Selleck Chemicals) at 33°C for up to 3 days. Disseminations of cells were monitored using an Olympus SZX10 fluorescent stereo microscope equipped with an Olympus DP71 camera.

Statistical tests

We compared data from Western blots, apoptosis assays, and zebrafish xenografts using nonparametric hypothesis tests, because these tests do not make any tacit assumption about the underlying distribution of the data. Kruskal-Wallis test (63) was used to compare data from two or more conditions (for example, control, drug treatment, and siRNA transfection), whereas Friedman’s test (64) was used to compare data that have multiple covariates (for example, time course measurements in control and treated samples, the two covariates being time points and treatment).

SUPPLEMENTARY MATERIALS

www.sciencesignaling.org/cgi/content/full/9/455/ra114/DC1

Note S1. Details of the BMRA algorithm.

Note S2. Benchmarking the BMRA algorithm.

Note S3. Sensitivity of the BMRA algorithm to noise and data incompleteness.

Note S4. Mechanistic modeling and parameter calibration using the BMM algorithm.

Note S5. Implementing BMRA to reconstruct the EGFR and IGF1R pathways.

Note S6. Mechanistic modeling of the EGFR and IGF1R pathways and the BMM algorithm.

Fig. S1. Interaction strengths of EGF and IGF1 network based on model simulation and experimental data.

Fig. S2. Effect of AKT inhibitor on IRS1 phosphorylation in four CRC cell lines.

Fig. S3. Effect of p70S6K knockdown on IRS1 phosphorylation in four CRC cell lines.

Fig. S4. Dynamics of ERK and AKT phosphorylation in HCT116 cells exposed to different concentrations of ligand.

Fig. S5. Effects of p70S6K knockdown and BIBX1382 treatment on serum-grown HCT116 cells.

Fig. S6. The amount of phosphorylated ERK in response to p70S6K knockdown in serum-grown HCT116 cells.

Fig. S7. Combined effect of p70S6K knockdown and BIBX1382 treatment on apoptosis of HCT116 cells.

Fig. S8. Combined effect of p70S6K knockdown and BIBX1382 treatment on apoptosis (cell death) of HKE3 cells.

Fig. S9. Combined effect of p70S6K knockdown and BIBX1382 treatment on apoptosis (cell death) of HT29 cells.

Code S1. MATLAB source codes for the BMRA, BMM algorithm, and the ODE models.

References (65–81)

REFERENCES AND NOTES

  1. ↵
    1. W. Kolch,
    2. M. Halasz,
    3. M. Granovskaya,
    4. B. N. Kholodenko
    , The dynamic control of signal transduction networks in cancer cells. Nat. Rev. Cancer 15, 515–527 (2015).
    OpenUrlCrossRefPubMed
  2. ↵
    1. F. Bertaux,
    2. S. Stoma,
    3. D. Drasdo,
    4. G. Batt
    , Modeling dynamics of cell-to-cell variability in TRAIL-induced apoptosis explains fractional killing and predicts reversible resistance. PLOS Comput. Biol. 10, e1003893 (2014).
    OpenUrl
  3. ↵
    1. N. Borisov,
    2. E. Aksamitiene,
    3. A. Kiyatkin,
    4. S. Legewie,
    5. J. Berkhout,
    6. T. Maiwald,
    7. N. P. Kaimachnikov,
    8. J. Timmer,
    9. J. B. Hoek,
    10. B. N. Kholodenko
    , Systems-level interactions between insulin–EGF networks amplify mitogenic signaling. Mol. Syst. Biol. 5, 256 (2009).
    OpenUrlAbstract/FREE Full Text
    1. C. Chen,
    2. W. T. Baumann,
    3. R. Clarke,
    4. J. J. Tyson
    , Modeling the estrogen receptor to growth factor receptor signaling switch in human breast cancer cells. FEBS Lett. 587, 3327–3334 (2013).
    OpenUrlCrossRefPubMedWeb of Science
    1. C. Chen,
    2. W. T. Baumann,
    3. J. Xing,
    4. L. Xu,
    5. R. Clarke,
    6. J. J. Tyson
    , Mathematical models of the transitions between endocrine therapy responsive and resistant states in breast cancer. J. R. Soc. Interface 11, 20140206 (2014).
    OpenUrlAbstract/FREE Full Text
  6. ↵
    1. D. Fey,
    2. D. R. Croucher,
    3. W. Kolch,
    4. B. N. Kholodenko
    , Crosstalk and signaling switches in mitogen-activated protein kinase cascades. Front. Physiol. 3, 355 (2012).
    OpenUrlCrossRefPubMed
    1. A. Gambin,
    2. A. Charzyńska,
    3. A. Ellert-Miklaszewska,
    4. M. Rybiński
    , Computational models of the JAK1/2-STAT1 signaling. JAKSTAT 2, e24672 (2013).
    OpenUrl
  8. ↵
    1. B. N. Kholodenko,
    2. O. V. Demin,
    3. G. Moehren,
    4. J. B. Hoek
    , Quantification of short term signaling by the epidermal growth factor receptor. J. Biol. Chem. 274, 30169–30181 (1999).
    OpenUrlAbstract/FREE Full Text
  9. ↵
    1. R. J. Orton,
    2. M. E. Adriaens,
    3. A. Gormand,
    4. O. E. Sturm,
    5. W. Kolch,
    6. D. R. Gilbert
    , Computational modelling of cancerous mutations in the EGFR/ERK signalling pathway. BMC Syst. Biol. 3, 100 (2009).
    OpenUrlCrossRefPubMed
  10. ↵
    1. B. B. Aldridge,
    2. J. M. Burke,
    3. D. A. Lauffenburger,
    4. P. K. Sorger
    , Physicochemical modelling of cell signalling pathways. Nat. Cell Biol. 8, 1195–1203 (2006).
    OpenUrlCrossRefPubMedWeb of Science
  11. ↵
    1. B. N. Kholodenko
    , Cell-signalling dynamics in time and space. Nat. Rev. Mol. Cell Biol. 7, 165–176 (2006).
    OpenUrlCrossRefPubMedWeb of Science
  12. ↵
    1. P. K. Kreeger,
    2. R. Mandhana,
    3. S. K. Alford,
    4. K. M. Haigis,
    5. D. A. Lauffenburger
    , RAS mutations affect tumor necrosis factor–induced apoptosis in colon carcinoma cells via ERK-modulatory negative and positive feedback circuits along with non-ERK pathway effects. Cancer Res. 69, 8191–8199 (2009).
    OpenUrlAbstract/FREE Full Text
    1. M. K. Morris,
    2. J. Saez-Rodriguez,
    3. D. C. Clarke,
    4. P. K. Sorger,
    5. D. A. Lauffenburger
    , Training signaling pathway maps to biochemical data with constrained fuzzy logic: Quantitative analysis of liver cell responses to inflammatory stimuli. PLOS Comput. Biol. 7, e1001099 (2011).
    OpenUrlCrossRefPubMed
    1. A. R. Tentner,
    2. M. J. Lee,
    3. G. J. Ostheimer,
    4. L. D. Samson,
    5. D. A. Lauffenburger,
    6. M. B. Yaffe
    , Combined experimental and computational analysis of DNA damage signaling reveals context-dependent roles for Erk in apoptosis and G1/S arrest after genotoxic stress. Mol. Syst. Biol. 8, 568 (2012).
    OpenUrlAbstract/FREE Full Text
  15. ↵
    1. S. Mukherjee,
    2. T. P. Speed
    , Network inference using informative priors. Proc. Natl. Acad. Sci. U.S.A. 105, 14313–14318 (2008).
    OpenUrlAbstract/FREE Full Text
  16. ↵
    1. K. Sachs,
    2. O. Perez,
    3. D. Pe’er,
    4. D. A. Lauffenburger,
    5. G. P. Nolan
    , Causal protein-signaling networks derived from multiparameter single-cell data. Science 308, 523–529 (2005).
    OpenUrlAbstract/FREE Full Text
  17. ↵
    1. Ö. Sahin,
    2. H. Fröhlich,
    3. C. Löbke,
    4. U. Korf,
    5. S. Burmester,
    6. M. Majety,
    7. J. Mattern,
    8. I. Schupp,
    9. C. Chaouiya,
    10. D. Thieffry,
    11. A. Poustka,
    12. S. Wiemann,
    13. T. Beissbarth,
    14. D. Arlt
    , Modeling ERBB receptor-regulated G1/S transition to find novel targets for de novo trastuzumab resistance. BMC Syst. Biol. 3, 1 (2009).
    OpenUrlCrossRefPubMed
  18. ↵
    1. T. Santra,
    2. W. Kolch,
    3. B. N. Kholodenko
    , Integrating Bayesian variable selection with modular response analysis to infer biochemical network topology. BMC Syst. Biol. 7, 57 (2013).
    OpenUrl
  19. ↵
    1. B. N. Kholodenko,
    2. A. Kiyatkin,
    3. F. J. Bruggeman,
    4. E. Sontag,
    5. H. V. Westerhoff,
    6. J. B. Hoek
    , Untangling the wires: A strategy to trace functional interactions in signaling and gene networks. Proc. Natl. Acad. Sci. U.S.A. 99, 12841–12846 (2002).
    OpenUrlAbstract/FREE Full Text
  20. ↵
    1. M. Andrec,
    2. B. N. Kholodenko,
    3. R. M. Levy,
    4. E. Sontag
    , Inference of signaling and gene regulatory networks by steady-state perturbation experiments: Structure and accuracy. J. Theor. Biol. 232, 427–441 (2005).
    OpenUrlCrossRefPubMedWeb of Science
    1. S. D. M. Santos,
    2. P. J. Verveer,
    3. P. I. H. Bastiaens
    , Growth factor-induced MAPK network topology shapes Erk response determining PC-12 cell fate. Nat. Cell Biol. 9, 324–330 (2007).
    OpenUrlCrossRefPubMedWeb of Science
  22. ↵
    1. B. Klinger,
    2. A. Sieber,
    3. R. Fritsche-Guenther,
    4. F. Witzel,
    5. L. Berry,
    6. D. Schumacher,
    7. Y. Yan,
    8. P. Durek,
    9. M. Merchant,
    10. R. Schäfer,
    11. C. Sers,
    12. N. Blüthgen
    , Network quantification of EGFR signaling unveils potential for targeted combination therapy. Mol. Syst. Biol. 9, 673 (2013).
    OpenUrlAbstract/FREE Full Text
  23. ↵
    1. A. A. Margolin,
    2. I. Nemenman,
    3. K. Basso,
    4. C. Wiggins,
    5. G. Stolovitzky,
    6. R. D. Favera,
    7. A. Califano
    , ARACNE: An algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context. BMC Bioinformatics 7 (Suppl. 1), S7 (2006).
    OpenUrlCrossRefPubMed
  24. ↵
    1. M. F. Ciaccio,
    2. J. P. Wagner,
    3. C.-P. Chuu,
    4. D. A. Lauffenburger,
    5. R. B. Jones
    , Systems analysis of EGF receptor signaling dynamics with microwestern arrays. Nat. Methods 7, 148–155 (2010).
    OpenUrlCrossRefPubMedWeb of Science
  25. ↵
    1. C. Siegenthaler,
    2. R. Gunawan
    , Assessment of network inference methods: How to cope with an underdetermined problem. PLOS ONE 9, e90481 (2014).
    OpenUrlCrossRefPubMed
  26. ↵
    1. M. F. Ciaccio,
    2. V. C. Chen,
    3. R. B. Jones,
    4. N. Bagheri
    , The DIONESUS algorithm provides scalable and accurate reconstruction of dynamic phosphoproteomic networks to reveal new drug targets. Integr. Biol. 7, 776–791 (2015).
    OpenUrl
    1. Y. Şenbabaoğlu,
    2. S. O. Sümer,
    3. F. Sánchez-Vega,
    4. D. Bemis,
    5. G. Ciriello,
    6. N. Schultz,
    7. C. Sander
    , A multi-method approach for proteomic network inference in 11 human cancers. PLOS Comput. Biol. 12, e1004765 (2016).
    OpenUrl
    1. T. Qin,
    2. L. C. Tsoi,
    3. K. J. Sims,
    4. X. Lu,
    5. W. J. Zheng
    , Signaling network prediction by the Ontology Fingerprint enhanced Bayesian network. BMC Syst. Biol. 6 (Suppl. 3), S3 (2012).
    OpenUrl
  29. ↵
    1. Z. Wang,
    2. W. Xu,
    3. F. A. San Lucas,
    4. Y. Liu
    , Incorporating prior knowledge into Gene Network Study. Bioinformatics 29, 2633–2640 (2013).
    OpenUrlAbstract/FREE Full Text
  30. ↵
    1. B. Asadi,
    2. M. R. Maurya,
    3. D. M. Tartakovsky,
    4. S. Subramaniam
    , Comparison of statistical and optimisation-based methods for data-driven network reconstruction of biochemical systems. IET Syst. Biol. 6, 155–163 (2012).
    OpenUrlPubMed
    1. C. Cosentino,
    2. W. Curatola,
    3. F. Montefusco,
    4. M. Bansal,
    5. D. di Bernardo,
    6. F. Amato
    , Linear matrix inequalities approach to reconstruction of biological networks. IET Syst. Biol. 1, 164–173 (2007).
    OpenUrlPubMed
  32. ↵
    1. F. Farhangmehr,
    2. M. R. Maurya,
    3. D. M. Tartakovsky,
    4. S. Subramaniam
    , Information theoretic approach to complex biological network reconstruction: Application to cytokine release in RAW 264.7 macrophages. BMC Syst. Biol. 8, 77 (2014).
    OpenUrl
  33. ↵
    1. T. Grote,
    2. D. R. Siwak,
    3. H. A. Fritsche,
    4. C. Joy,
    5. G. B. Mills,
    6. D. Simeone,
    7. D. C. Whitcomb,
    8. C. D. Logsdon
    , Validation of reverse phase protein array for practical screening of potential biomarkers in serum and plasma: Accurate detection of CA19-9 levels in pancreatic cancer. Proteomics 8, 3051–3060 (2008).
    OpenUrlCrossRefPubMedWeb of Science
  34. ↵
    1. A. M. Krasinskas
    , EGFR signaling in colorectal carcinoma. Patholog. Res. Int. 2011, 932932 (2011).
    OpenUrlCrossRefPubMed
  35. ↵
    1. W. De Roock,
    2. B. Claes,
    3. D. Bernasconi,
    4. J. De Schutter,
    5. B. Biesmans,
    6. G. Fountzilas,
    7. K. T. Kalogeras,
    8. V. Kotoula,
    9. D. Papamichael,
    10. P. Laurent-Puig,
    11. F. Penault-Llorca,
    12. P. Rougier,
    13. B. Vincenzi,
    14. D. Santini,
    15. G. Tonini,
    16. F. Cappuzzo,
    17. M. Frattini,
    18. F. Molinari,
    19. P. Saletti,
    20. S. De Dosso,
    21. M. Martini,
    22. A. Bardelli,
    23. S. Siena,
    24. A. Sartore-Bianchi,
    25. J. Tabernero,
    26. T. Macarulla,
    27. F. Di Fiore,
    28. A. O. Gangloff,
    29. F. Ciardiello,
    30. P. Pfeiffer,
    31. C. Qvortrup,
    32. T. P. Hansen,
    33. E. Van Cutsem,
    34. H. Piessevaux,
    35. D. Lambrechts,
    36. M. Delorenzi,
    37. S. Tejpar
    , Effects of KRAS, BRAF, NRAS, and PIK3CA mutations on the efficacy of cetuximab plus chemotherapy in chemotherapy-refractory metastatic colorectal cancer: A retrospective consortium analysis. Lancet Oncol. 11, 753–762 (2010).
    OpenUrlCrossRefPubMedWeb of Science
  36. ↵
    1. B. O. Van Emburgh,
    2. A. Sartore-Bianchi,
    3. F. Di Nicolantonio,
    4. S. Siena,
    5. A. Bardelli
    , Acquired resistance to EGFR-targeted therapies in colorectal cancer. Mol. Oncol. 8, 1084–1094 (2014).
    OpenUrlCrossRefPubMed
  37. ↵
    1. I. Singh,
    2. H. Amin,
    3. B. Rah,
    4. A. Goswami
    , Targeting EGFR and IGF 1R: A promising combination therapy for metastatic cancer. Front. Biosci. (Schol. Ed.) 5, 231–246 (2013).
    OpenUrl
  38. ↵
    1. T. Troiani,
    2. S. Napolitano,
    3. D. Vitagliano,
    4. F. Morgillo,
    5. A. Capasso,
    6. V. Sforza,
    7. A. Nappi,
    8. D. Ciardiello,
    9. F. Ciardiello,
    10. E. Martinelli
    , Primary and acquired resistance of colorectal cancer cells to anti-EGFR antibodies converge on MEK/ERK pathway activation and can be overcome by combined MEK/EGFR inhibition. Clin. Cancer Res. 25, 3775–3786 (2014).
    OpenUrl
  39. ↵
    1. E. Rozengurt,
    2. H. P. Soares,
    3. J. Sinnet-Smith
    , Suppression of feedback loops mediated by PI3K/mTOR induces multiple overactivation of compensatory pathways: An unintended consequence leading to drug resistance. Mol. Cancer Ther. 13, 2477–2488 (2014).
    OpenUrlAbstract/FREE Full Text
  40. ↵
    1. S. Shirasawa,
    2. M. Furuse,
    3. N. Yokoyama,
    4. T. Sasazuki
    , Altered growth of human colon cancer cell lines disrupted at activated Ki-ras. Science 260, 85–88 (1993).
    OpenUrlAbstract/FREE Full Text
  41. ↵
    1. E. Aksamitiene,
    2. A. Kiyatkin,
    3. B. N. Kholodenko
    , Cross-talk between mitogenic Ras/MAPK and survival PI3K/Akt pathways: A fine balance. Biochem. Soc. Trans. 40, 139–146 (2012).
    OpenUrlAbstract/FREE Full Text
  42. ↵
    1. M. A. Tsyganov,
    2. W. Kolch,
    3. B. N. Kholodenko
    , The topology design principles that determine the spatiotemporal dynamics of G-protein cascades. Mol. Biosyst. 8, 730–743 (2012).
    OpenUrlCrossRefPubMed
  43. ↵
    1. S. M. Johnson,
    2. P. Gulhati,
    3. I. Arrieta,
    4. X. Wang,
    5. T. Uchida,
    6. T. Gao,
    7. B. M. Evers
    , Curcumin inhibits proliferation of colorectal carcinoma by modulating Akt/mTOR signaling. Anticancer Res. 29, 3185–3190 (2009).
    OpenUrlAbstract/FREE Full Text
  44. ↵
    1. R. S. Sawhney,
    2. M. M. Cookson,
    3. B. Sharma,
    4. J. Hauser,
    5. M. G. Brattain
    , Autocrine transforming growth factor α regulates cell adhesion by multiple signaling via specific phosphorylation sites of p70S6 kinase in colon cancer cells. J. Biol. Chem. 279, 47379–47390 (2004).
    OpenUrlAbstract/FREE Full Text
  45. ↵
    1. E. F. Dunn,
    2. M. Iida,
    3. R. A. Myers,
    4. D. A. Campbell,
    5. K. A. Hintz,
    6. E. A. Armstrong,
    7. C. Li,
    8. D. L. Wheeler
    , Dasatinib sensitizes KRAS mutant colorectal tumors to cetuximab. Oncogene 30, 561–574 (2011).
    OpenUrlCrossRefPubMed
    1. E. Aksamitiene,
    2. B. N. Kholodenko,
    3. W. Kolch,
    4. J. B. Hoek,
    5. A. Kiyatkin
    , PI3K/Akt-sensitive MEK-independent compensatory circuit of ERK activation in ER-positive PI3K-mutant T47D breast cancer cells. Cell. Signal. 22, 1369–1378 (2010).
    OpenUrlCrossRefPubMed
    1. Q. Ye,
    2. W. Cai,
    3. Y. Zheng,
    4. B. M. Evers,
    5. Q.-B. She
    , ERK and AKT signaling cooperate to translationally regulate survivin expression for metastatic progression of colorectal cancer. Oncogene 33, 1828–1839 (2014).
    OpenUrlCrossRefPubMed
    1. Q. G. Jiang,
    2. T. Y. Li,
    3. D. N. Liu,
    4. H. T. Zhang
    , PI3K/Akt pathway involving into apoptosis and invasion in human colon cancer cells LoVo. Mol. Biol. Rep. 41, 3359–3367 (2014).
    OpenUrl
    1. J.-J. Liu,
    2. R.-D. Duan
    , LY294002 enhances boswellic acid-induced apoptosis in colon cancer cells. Anticancer Res 29, 2987–2991 (2009).
    OpenUrlAbstract/FREE Full Text
  50. ↵
    1. S. Ponnurangam,
    2. D. Standing,
    3. P. Rangarajan,
    4. D. Subramaniam
    , Tandutinib inhibits the Akt/mTOR signaling pathway to inhibit colon cancer growth. Mol. Cancer Ther. 12, 598–609 (2013).
    OpenUrlAbstract/FREE Full Text
  51. ↵
    1. M. Girolami,
    2. B. Calderhead
    , Riemann manifold Langevin and Hamiltonian Monte Carlo methods. J. R. Stat. Soc. Series B Stat. Methodology 73, 123–214 (2011).
    OpenUrl
  52. ↵
    1. D. Chaturvedi,
    2. X. Gao,
    3. M. S. Cohen,
    4. J. Taunton,
    5. T. B. Patel
    , Rapamycin induces transactivation of the EGFR and increases cell survival. Oncogene 28, 1187–1196 (2009).
    OpenUrlCrossRefPubMedWeb of Science
  53. ↵
    1. E. Buck,
    2. A. Eyzaguirre,
    3. M. Rosenfeld-Franklin,
    4. S. Thomson,
    5. M. Mulvihill,
    6. S. Barr,
    7. E. Brown,
    8. M. O’Connor,
    9. Y. Yao,
    10. J. Pachter,
    11. M. Miglarese,
    12. D. Epstein,
    13. K. K. Iwata,
    14. J. D. Haley,
    15. N. W. Gibson,
    16. Q.-S. Ji
    , Feedback mechanisms promote cooperativity for small molecule inhibitors of epidermal and insulin-like growth factor receptors. Cancer Res. 68, 8322–8332 (2008).
    OpenUrlAbstract/FREE Full Text
  54. ↵
    1. K. Schmid,
    2. Z. Bago-Horvath,
    3. W. Berger,
    4. A. Haitel,
    5. D. Cejka,
    6. J. Werzowa,
    7. M. Filipits,
    8. B. Herberger,
    9. H. Hayden,
    10. W. Sieghart
    , Dual inhibition of EGFR and mTOR pathways in small cell lung cancer. Br. J. Cancer 103, 622–628 (2010).
    OpenUrlCrossRefPubMedWeb of Science
  55. ↵
    1. N. E. Bhola,
    2. S. M. Thomas,
    3. M. Freilino,
    4. S. Joyce,
    5. A. Sahu,
    6. J. Maxwell,
    7. A. Argiris,
    8. R. Seethala,
    9. J. R. Grandis
    , Targeting GPCR-mediated p70S6K activity may improve head and neck cancer response to cetuximab. Clin. Cancer Res. 17, 4996–5004 (2011).
    OpenUrlAbstract/FREE Full Text
  56. ↵
    1. M. Gupta,
    2. J. G. Ibrahim
    , An information matrix prior for Bayesian analysis in generalized linear models with high dimensional data. Stat. Sin. 19, 1641–1663 (2009).
    OpenUrlPubMed
  57. ↵
    1. T. Bayes,
    2. R. Price
    , An essay towards solving a problem in the doctrine of chances. By the late Rev. Mr. Bayes, F. R. S. communicated by Mr. Price, in a letter to John Canton, A. M. F. R. S. Philos. Trans. 53, 370–418 (1763).
    OpenUrlFREE Full Text
  58. ↵
    1. T. Santra
    , A Bayesian framework that integrates heterogeneous data for inferring gene regulatory networks. Front. Bioeng. Biotechnol. 2, 13 (2014).
    OpenUrl
  59. ↵
    1. W. K. Hastings
    , Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57, 97–109 (1970).
    OpenUrlAbstract/FREE Full Text
  60. ↵
    1. H. Haario,
    2. E. Saksman,
    3. J. Tamminen
    , An adaptive Metropolis algorithm. Bernoulli 7, 223–242 (2001).
    OpenUrlCrossRefWeb of Science
  61. ↵
    C. J. Geyer, in Proceedings of the 23rd Symposium on the Interface, E. M. Keramidas, S. M. Kaufman, Eds. (American Statistical Association, 1991), p. 156.
  62. ↵
    1. M. Halasz,
    2. B. Polgar,
    3. G. Berta,
    4. L. Czimbalek,
    5. J. Szekeres-Bartho
    , Progesterone-induced blocking factor differentially regulates trophoblast and tumor invasion by altering matrix metalloproteinase activity. Cell. Mol. Life Sci. 70, 4617–4630 (2013).
    OpenUrl
  63. ↵
    1. W. H. Kruskal,
    2. W. A. Wallis
    , Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 47, 583–621 (1952).
    OpenUrlCrossRefWeb of Science
  64. ↵
    1. M. Friedman
    , The use of ranks to avoid the assumption of normality implicit in the analysis of variance. J. Am. Stat. Assoc. 32, 675–701 (1937).
    OpenUrlCrossRefWeb of Science
  65. ↵
    1. D. Marbach,
    2. R. J. Prill,
    3. T. Schaffter,
    4. C. Mattiussi,
    5. D. Floreano,
    6. G. Stolovitzky
    , Revealing strengths and weaknesses of methods for gene network inference. Proc. Natl. Acad. Sci. U.S.A. 107, 6286–6291 (2010).
    OpenUrlAbstract/FREE Full Text
    1. D. Marbach,
    2. T. Schaffter,
    3. C. Mattiussi,
    4. D. Floreano
    , Generating realistic in silico gene networks for performance assessment of reverse engineering methods. J. Comput. Biol. 16, 229–239 (2009).
    OpenUrlCrossRefPubMedWeb of Science
    1. D. M. W. Powers
    , Evaluation: From precision, recall and f-measure to roc., informedness, markedness & correlation. J. Mach. Learn. Technol. 2, 37–63 (2011).
    OpenUrl
    1. F. Montefusco,
    2. C. Cosentino,
    3. F. Amato
    , CORE-Net: Exploiting prior knowledge and preferential attachment to infer biological interaction networks. IET Syst. Biol. 4, 296–310 (2010).
    OpenUrlPubMed
    1. K. Hiratani,
    2. T. Haruta,
    3. A. Tani,
    4. J. Kawahara,
    5. I. Usui,
    6. M. Kobayashi
    , Roles of mTOR and JNK in serine phosphorylation, translocation, and degradation of IRS-1. Biochem. Biophys. Res. Commun. 335, 836–842 (2005).
    OpenUrlCrossRefPubMedWeb of Science
    1. C. L. Mamay,
    2. A. M. Mingo-Sion,
    3. D. M. Wolf,
    4. M. D. Molina,
    5. C. L. Van Den Berg
    , An inhibitory function for JNK in the regulation of IGF-I signaling in breast cancer. Oncogene 22, 602–614 (2003).
    OpenUrlCrossRefPubMedWeb of Science
    1. R. Roskoski Jr
    , MEK1/2 dual-specificity protein kinases: Structure and regulation. Biochem. Biophys. Res. Commun. 417, 5–10 (2012).
    OpenUrlCrossRefPubMed
    1. K. Moelling,
    2. K. Schad,
    3. M. Bosse,
    4. S. Zimmermann,
    5. M. Schweneker
    , Regulation of Raf-Akt Cross-talk. J. Biol. Chem. 277, 31099–31106 (2002).
    OpenUrlAbstract/FREE Full Text
    1. Q. Ding,
    2. W. Xia,
    3. J.-C. Liu,
    4. J.-Y. Yang,
    5. D.-F. Lee,
    6. J. Xia,
    7. G. Bartholomeusz,
    8. Y. Li,
    9. Y. Pan,
    10. Z. Li,
    11. R. C. Bargou,
    12. J. Qin,
    13. C.-C. Lai,
    14. F.-J. Tsai,
    15. C.-H. Tsai,
    16. M.-C. Hung
    , Erk associates with and primes GSK-3β for its inactivation resulting in upregulation of β-catenin. Mol. Cell 19, 159–170 (2005).
    OpenUrlCrossRefPubMedWeb of Science
    1. Y. Shi,
    2. J.-h. Hsu,
    3. L. Hu,
    4. J. Gera,
    5. A. Lichtenstein
    , Signal pathways involved in activation of p70S6K and phosphorylation of 4E-BP1 following exposure of multiple myeloma tumor cells to interleukin-6. J. Biol. Chem. 277, 15712–15720 (2002).
    OpenUrlAbstract/FREE Full Text
    1. B.-C. Chen,
    2. W.-W. Lin
    , PKC- and ERK-dependent activation of IκB kinase by lipopolysaccharide in macrophages: Enhancement by P2Y receptor-mediated CaMK activation. Br. J. Pharmacol. 134, 1055–1065 (2001).
    OpenUrlCrossRefPubMedWeb of Science
    1. X. Fang,
    2. S. X. Yu,
    3. Y. Lu,
    4. R. C. Bast Jr,
    5. J. R. Woodgett,
    6. G. B. Mills
    , Phosphorylation and inactivation of glycogen synthase kinase 3 by protein kinase A. Proc. Natl. Acad. Sci. U.S.A. 97, 11960–11965 (2000).
    OpenUrlAbstract/FREE Full Text
    1. D. A. E. Cross,
    2. D. R. Alessi,
    3. P. Cohen,
    4. M. Andjelkovich,
    5. B. A. Hemmings
    , Inhibition of glycogen synthase kinase-3 by insulin mediated by protein kinase B. Nature 378, 785–789 (1995).
    OpenUrlCrossRefPubMedWeb of Science
    1. D. Bai,
    2. L. Ueno,
    3. P. K. Vogt
    , Akt-mediated regulation of NFκB and the essentialness of NFκB for the oncogenicity of PI3K and Akt. Int. J. Cancer 125, 2863–2870 (2009).
    OpenUrlCrossRefPubMedWeb of Science
    1. H. P. Soares,
    2. Y. Ni,
    3. K. Kisfalvi,
    4. J. Sinnett-Smith,
    5. E. Rozengurt
    , Different patterns of Akt and ERK feedback activation in response to rapamycin, active-site mTOR inhibitors and metformin in pancreatic cancer cells. PLOS ONE 8, e57289 (2013).
    OpenUrlCrossRefPubMed
    1. S. Y. Yoon,
    2. Y. J. Lee,
    3. J. H. Seo,
    4. H. J. Sung,
    5. K. H. Park,
    6. I. K. Choi,
    7. S. J. Kim,
    8. S. C. Oh,
    9. C. W. Choi,
    10. B. S. Kim,
    11. S. W. Shin,
    12. Y. H. Kim,
    13. J. S. Kim
    , uPAR expression under hypoxic conditions depends on iNOS modulated ERK phosphorylation in the MDA-MB-231 breast carcinoma cell line. Cell Res. 16, 75–81 (2006).
    OpenUrlCrossRefPubMedWeb of Science
  81. ↵
    1. B. Kholodenko,
    2. M. B. Yaffe,
    3. W. Kolch
    , Computational approaches for analyzing information flow in biological networks. Sci. Signal. 5, re1 (2012).
    OpenUrlAbstract/FREE Full Text
Acknowledgments: We thank D. Fey for stimulating discussions on several aspects of the computational work. We also thank A. Blanco (Conway Core Technologies, University College Dublin) for his help with flow cytometry and B. Kennedy (Zebrafish Facility, University College Dublin) for his useful advice on microinjections. We thank S. Shirasawa for providing us the HCT116 and HKE3 CRC cell lines. Funding: This work was supported by Science Foundation Ireland grant 06/CE/B1129, European Union FP7 grant 259348-2 “ASSET,” and Irish Cancer Society Collaborative Cancer Research Centre BREAST-PREDICT grant CCRC13GAL. B.N.K. acknowledges support from the EU FP7 SynSignal (Grant agreement number 613879) and EU H2020 SmartNanoTox (Grant agreement number 686098). Author contributions: M.H. designed and performed the biological experiments and wrote the paper; B.N.K. designed the study and wrote the paper; W.K. designed the study and wrote the paper; T.S. designed the study, designed and implemented the computational analysis, and wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The mathematical model used in this study is available as part of the Supplementary Materials.
View Abstract

Article Tools

My saved folders

Stay Connected to Science Signaling

Related Content

Similar Articles in:

Citing Articles in:

Navigate This Article