## 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).

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.

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 Ser^{473}, extracellular signal–regulated kinase 2 (ERK2) Thr^{185}/Tyr^{187}, MAPK kinase 1 (MEK1) Ser^{217}/Ser^{221}, p70S6K Thr^{421}/Ser^{424}, IGF1R Tyr^{1131}, glycogen synthase kinase 3 α/β (GSK3α/β) Ser^{21}/Ser^{9}, IκBα Ser^{32}/Ser^{36}, and insulin receptor substrate 1 (IRS1) Ser^{636}/Ser^{639}] 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).

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 (Ser^{636} and Ser^{639}). 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 IRS1^{S636/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 IRS1^{S636/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 IRS1^{S636/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 IRS1^{S636/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.

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).

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 (*k*_{ras1} = 1.35 ± 0.36 min^{− 1} and *k*_{ras2} = 3.02 ± 1.5 min^{−1}; see Fig. 5C) than when it is inactivated [*V*_{Mras} = 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*).

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.

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 log_{2} 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).

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 (*n*_{i}, *i* = 1 … *N*) biochemical network. MRA (*19*) showed that following a perturbation (*p*_{k}), the fractional change *x*_{i}^{SS}) of a network node (*n*_{i}) is linearly related to that (*R*_{jk}, *j* ≠ *i*) of the other nodes (*n*_{j}, *j* ≠ *i*) through the corresponding interaction strengths (*r*_{ij}, *j* ≠ *i*).

If node *n*_{j}, *j* ≠ *i* directly influences node *n*_{i}, then *r*_{ij} ≠ 0; otherwise, *r*_{ij} = 0. We modified this equation by introducing two more variables, a binary variable *A*_{ij} in Eq. 1 to indicate whether *n*_{j} directly influences *n*_{i} (*A*_{ij} = 1) or not (*A*_{ij} = 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.

We used Bayesian statistics to infer the probability distributions of *A*_{ij}, *i* ≠ *j* and *r*_{ij}, *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 (*n*_{i}) and its regulators. The probability distributions of the network descriptors (**A _{i}** = {

*A*

_{ij},

*j*≠

*i*},

**r**= {

_{i}*r*

_{ij},

*j*≠

*i*}) of each of these subnetworks were inferred separately. We first formulated prior probability distributions of

**A**and

_{i}**r**based on existing knowledge of the network topology. The binary vector

_{i}*n*

_{i}. The prior distribution (

**P**(

**A**)) of

_{i}**A**was then formulated as follows:

_{i}**A**and

_{i}*P*(

**A**) 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

_{i}**r**is dependent on

_{i}**A**and is denoted by

_{i}*P*(

**r**|

_{i}**A**). In the absence of a direct interaction from

_{i}*n*

_{j}to

*n*

_{i}(that is,

*A*

_{ij}= 0), the corresponding interaction strength (

*r*

_{ij}) was assumed to have 0 value with probability 1, whereas the interaction strengths representing direct interactions (

*A*

_{ij}= 1,

*j*≠

*i*) were assumed to have Gaussian priors (note S1) (

*18*,

*56*). The above prior distributions (

*P*(

**A**),

_{i}*P*(

**r**|

_{i}**A**)) were then used to calculate the likelihood (

_{i}*P*(

**R**|

**r**,

_{i}**A**)) of the global changes (

_{i}**R**) given the network descriptors (

**A**,

_{i}**r**) (note S1). Subsequently, we applied Bayes’ rule (

_{i}*57*) to calculate the joint posterior distribution (

*P*(

**A**,

_{i}**r**|

_{i}**)) for the network descriptors (**

*R***r**,

_{i}**A**). Using the chain rule, the above joint posterior

_{i}*P*(

**A**,

_{i}**r**|

_{i}**R**)) is the product of the marginal posterior (

*P*(

**A**|

_{i}**R**)) of the binary indicators (

**A**) and the conditional posterior (

_{i}*P*(

**r**|

_{i}**A**,

_{i}**R**)) of the interactions strengths (

**r**), that is,

_{i}*P*(

**A**,

_{i}**r**|

_{i}**R**) =

*P*(

**r**|

_{i}**A**,

_{i}**R**) ×

*P*(

**A**|

_{i}**R**). The marginal posterior (

*P*(

**A**|

_{i}**R**)) is proportional to the product of a multivariate Student’s

*t*(MVSt) and the prior distribution, whereas the conditional posterior (

*P*(

**r**|

_{i}**A**,

_{i}**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*(

**A**,

_{i}**r**|

_{i}**R**). The MCMC algorithm starts with a random network, and each iteration (

*t*) randomly proposes a new subnetwork (

*15*,

*58*). The proposed subnetwork was then accepted (

*59*). Subsequently, we drew a sample (

^{5}) of iterations, and after an initial burn-in period (5 × 10

^{4}iterations), the accepted network topologies and interaction strengths were recorded. The mean of an element

*A*

_{ij}of

**A**over the post–burn-in samples was then used as the posterior probability of the interaction from node

_{i}*n*

_{j}to

*n*

_{i}, that is,

*r*

_{ij}) across the post–burn-in samples were used as the posterior mean, variance, and covariance, respectively, of the strengths of the interaction from node

*n*

_{j}to

*n*

_{i}. 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 (

*M*

_{o}) is defined as the area under the concentration versus time curve (

*x*

_{o}(

*t*)) of the output node (

*n*

_{o}) of the model, that is,

*n*

_{o1},

*n*

_{o2}…), the sum of the areas of their concentration curves (

*x*

_{o1}(

*t*),

*x*

_{o2}(

*t*), …) is used as the model output (

_{k}) in a parameter (θ

_{k}) causes the model output to change by Δ

*M*

_{o}, then the sensitivity (

*S*

_{k}) of the model to this parameter is given by

*S*

_{k}) were estimated for each model parameter (θ

_{k}) for a large number (10

^{5}) 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*(

**Θ**))(note S6) to the sensitive parameters (

_{s}**Θ**) and formulated a likelihood function (

_{s}**r**represents the simulated interaction strengths in condition

_{l}*l*(for example, ligand stimulation);

**μ**and

_{rl}**Σ**are their means and covariance matrices as estimated by the BMRA algorithm from experimental data. We then used Bayes’ rule (

_{rl}*57*) to calculate the posterior distribution [

*P*(

**Θ**|

_{s}**μ**,

_{rl}**Σ**,

_{rl}*l*= 1, 2 …)] of the sensitive parameters (

**Θ**):

_{s}*t*) proposes a new set of values (

*59*). This process was repeated many times (2.5 × 10

^{5}), and the samples drawn after a burn-in (10

^{5}) period were used to estimate the posterior distribution of

**Θ**The variances

_{s}.*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

**Θ**A MATLAB implementation of the above calibration algorithm can be found in the Supplementary Materials (code S1).

_{s}.### 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% CO_{2}.

### 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 × 10^{5} 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 (Ser^{473}), AKT, phospho-p44/42 MAPK (Thr^{202}/Tyr^{204}), p44/42 MAPK, p70 S6 kinase, phospho-p70 S6 kinase (Thr^{389}), IRS1, and phospho-IRS1 (Ser^{636}/Ser^{639}) 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 × 10^{5} 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-G_{1} 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 AND NOTES

**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.