## Abstract

Genomic features are used as biomarkers of sensitivity to kinase inhibitors used widely to treat human cancer, but effective patient stratification based on these principles remains limited in impact. Insofar as kinase inhibitors interfere with signaling dynamics, and, in turn, signaling dynamics affects inhibitor responses, we investigated associations in this study between cell-specific dynamic signaling pathways and drug sensitivity. Specifically, we measured 14 phosphoproteins under 43 different perturbed conditions (combinations of 5 stimuli and 7 inhibitors) in 14 colorectal cancer cell lines, building cell line–specific dynamic logic models of underlying signaling networks. Model parameters representing pathway dynamics were used as features to predict sensitivity to a panel of 27 drugs. Specific parameters of signaling dynamics correlated strongly with drug sensitivity for 14 of the drugs, 9 of which had no genomic biomarker. Following one of these associations, we validated a drug combination predicted to overcome resistance to MEK inhibitors by coblockade of GSK3, which was not found based on associations with genomic data. These results suggest that to better understand the cancer resistance and move toward personalized medicine, it is essential to consider signaling network dynamics that cannot be inferred from static genotypes. *Cancer Res; 77(12); 3364–75. ©2017 AACR*.

## Introduction

Patient response to anticancer therapies is extremely variable and understanding the reasons for this variability is a major challenge in cancer research. One approach to address this problem is to identify biomarkers that correlate with therapy response. However, except for a few examples, no efficient biomarkers are available (1). The problem of finding biomarkers has become even more important with the advent of targeted cancer therapies that are designed to affect specific molecular changes that drive the cancer, providing a larger number of treatment options. However, biomarkers that can be used to stratify patients for targeted drugs also remain largely elusive (2).

The most common approaches for patient stratification are currently based on genomic biomarkers, typically consisting of either copy number alteration or mutation of specific genes. The popularity of those biomarkers has been favored by the advancements of sequencing technologies and their subsequent decrease in cost. While for some drugs they have shown strong potential (3–5), for many other cases no efficient genomic marker for patient stratification exists, and for those where markers exist, they have rather low power due to the complex nature of cancer. Furthermore, their actual clinical significance for precision oncology is questionable (1).

Many targeted therapies inhibit specific signaling molecules, mostly kinases, and elicit a response that can affect and rewire large parts of the signaling network. Therefore, investigating the dynamics of signaling pathways can help to unveil new therapeutic strategies and biomarkers (6, 7). Understanding how drugs affect the signaling network as a whole, is particularly important as signaling dynamic differs not only between tumors of different tissue, but also between tumors of the same tissue. In fact, for different cancer cell types, the dynamics of signaling networks have been shown to differ strongly between different cells of the same tissue of origin (8–11).

Furthermore, models parameterized and experimentally calibrated using a specific cell line have been used to define pathway dynamic biomarkers of therapeutic outcome, such as drug sensitivity (12) or patient survival (13), based on model simulations in breast cancer and neuroblastoma, respectively.

The understanding of mechanisms of cellular response (network structure and dynamic behavior) can also be useful to tackle the development of compensatory signaling mechanisms of drug resistance (14), which is a recurrent problem for targeted therapies and is challenging to predict. Signal transduction is carried out by complex dynamic networks, and mutations in oncogenes or tumor suppressors can drastically change the properties of the network. This complex wiring confers robustness to the cells that often allows them to escape single-agent–targeted treatments (15, 16). On the basis of this idea, modeling of signaling pathways has been recently used to suggest combinatorial targeted therapies to effectively block multiple molecular pathways (8, 17–20): by integration of prior pathway knowledge with experimental observations, cell line–specific models were built, which could be used to simulate and thus prioritize possible combinatorial perturbation experiments.

In this article, we investigated to what extent dynamic interactions between different signaling pathways play a role in characterizing the specific cellular responses to drugs and how this can be used to suggest targeted combinatorial therapies. Furthermore, we assessed how these dynamic features of signaling fare against static genomic traits as markers of drug response. We did so by characterizing cell type–specific models for a panel of 14 different colorectal cancer cell lines that are integrated with a large-scale drug screening (21). We found associations with model parameters for 14 of the 27 drugs targeting our pathways of interests, for 9 of which there were no genomic biomarkers. These associations were used to define pathway dynamic biomarkers and to propose combinatorial therapies.

## Materials and Methods

### Cell lines

All cell lines used in this study (CAR1, CCK81, COLO-320-HSR, DIFI, HCT116, HT115, HT29, SK-CO-1, SNU-C5, SNU-C2B, SW1116, SW1463, SW620, and SW837) are mycoplasma negative. To exclude cross-contaminated or synonymous lines, a panel of 92 SNPs was profiled for each cell line (Sequenom) and a pairwise comparison score calculated for in-house identity checking across a panel of 1,000 cell lines (21). We confirmed the identity of cancer cell line set against those provided by the repositories, where possible, using a panel of 16 short tandem repeats (STR; AmpFLSTR Identifiler KIT, Applied Biosystems), which includes the 9 currently used by most of the cell line repositories (ATCC, Riken, JCRB, and DSMZ). STR or SNP data are available at COSMIC database (http://cancer.sanger.ac.uk/cell_lines). All cell lines are sourced from collaborators or repositories; they have been used for COSMIC cell line project and the Genomics of Drug Sensitivity in Cancer (GDSC) project (21) and have been stored in the laboratory as cryopreserved aliquots in liquid N_{2} for 7 years. A single cryovial is thawed for use and propagated for a maximum of 3 months before being discarded.

### Experimental setup of perturbation screen

Cell lines were cultured in either RPMI or DMEM/F12 medium with 10% FBS and were incubated at 37°C and 5% CO_{2}. Before perturbation commenced, cells were starved overnight in serum-free medium. At 90 minutes before lysis, the cells were treated with inhibitors (or solvent control DMSO) and at 30 minutes before lysis cells were stimulated with ligand (or solvent control H_{2}O). This procedure was conducted for all single inhibitor- and combined inhibitor–ligand perturbations. The following inhibitors were used: MEKi AZD6244 (4 μmol/L), PI3Ki AZD6482 (10 μmol/L), mTORi AZD8055 (2 μmol/L), TBK1i BX-795 (10 μmol/L), IKKi BMS-345541 (25 μmol/L), BRAFi PLX4720 (5 μmol/L), and TAK1i 5Z-7-Oxozeaenol (5 μmol/L). Concentrations were selected to inhibit target while minimizing off-target activity. We used the following ligands: EGF (25 ng/mL), HGF (50 ng/mL), IGF1 (10 ng/mL), TGFβ (5 ng/mL), and TNFα (10 ng/mL). After treatment and incubation, lysates were collected and analyzed with the Bio-Plex Protein Array system (Bio-Rad), as described in ref. 8, using magnetic beads specific for AKT^{S473}, c-Jun^{S63}, EGFR^{Y1068}, ERK1/2^{T202,Y204/T185,Y187}, GSK3A/B^{S21/S9}, IkBa^{S32,S36}, JNK^{T183,Y185}, MEK1^{S217,S221}, mTOR^{S2448}, p38^{T180,Y182}, PI3K^{Y458}, RPS6^{S235,S236}, p90RSK^{S380}, and SMAD2^{S465,S467}. The beads and detection antibodies were diluted 1:3. For data acquisition, the Bio-Plex Manager software and R package lxb was used.

### Data preprocessing

For each cell line, data were processed separately for each of the 14 measured phosphoproteins. The value of the control (unperturbed condition) was estimated as the median value of the 4 replicates and log2 fold changes with respect to the control were computed for each of the 42 perturbed conditions (data in Supplementary Table S1). As required by our logic formalism, data were linearly scaled between 0 and 1, with 0.5 corresponding to the basal (control) condition.

### Definition of logic ordinary differential equations

The logic ordinary differential equation formalism (22) is based on ordinary differential equations derived from logic models using a continuous update function *B _{i}* for each species (i.e., node in the network)

*x*, which can assume continuous values {0,1}. Differential equation for species

_{i}*i*is defined as follows:

Where *τ _{i}* is the life-time of species

*i*:

*τ*= 0 corresponds to a node not responsive to activation by upstream nodes and higher values of

_{i}*τ*represent faster response, interpretable as more functional species;

_{i}*x*are the

_{1,i}, x_{2,i}, …, x_{N,i}*N*regulators of

*x*and each regulation is defined by a sigmoidal transfer function

_{i}*f(x*, described in Supplementary Fig. S1, where parameter

_{i,j})*k*characterizes the strength of the regulation of species

_{j,i}*i*dependent on species

*j*(i.e., edge

*j → i*):

*k*= 0 means no regulation, while higher values represent stronger interactions.

_{j,i}Optimization of parameters *τ _{i}* and

*k*was performed using a modified version of the R package CNORode (add-on to CellNOptR; ref. 23) and the optimization toolbox MEIGO (24). To improve parameter estimates, we redefined the objective function (i.e., sum of the squared difference between model predictions and measured values,

_{j,i}*Q*) introducing an L1 penalty term on the parameters (L1 norm is the sum of the absolute values of the parameters, P). The balance between prioritizing good fit (low

_{LS}*Q*) or sparse model (low

_{LS}*P*) is regulated by a tunable parameter

*λ,*with new objective function:

In addition, the objective function includes a penalty to promote that the simulated values at *t* = 30 minutes (corresponding to the measured time point) are at steady state to match the experimental assumptions. The regularization was performed sequentially first on *τ _{i}* (to remove unresponsive nodes), and then, fixing

*τ*, on

_{i}*k*(to induce edge sparsity). At each step, the term

_{j,i}*λ*was tuned to guarantee a good compromise between goodness of fit and model sparsity (Supplementary Figs. S2 and S3). Finally, bootstrap distribution was obtained for parameters

*k*by repeating the optimization on resampled experimental data with replacement 150 times. More details are in Supplementary Methods, and code for modified version of the CNORode package at https://github.com/saezlab/CNORode2017.

_{j,i}### Elastic net

To select only robust features, Elastic net (25) training was iterated *M* times (*M* is the total number of observations), each time leaving out one observation (sampling without replacement). For each iteration, leave-one-out cross-validation was applied to tune Elastic Net hyperparameters (using *caret* R package). Median values were computed for the weights of the features across the *M* iterations, to select only features appearing in at least half of the iterations.

### Statistical analysis

ANOVA was used to define genomic markers of drug sensitivity using the GDSC tools (http://gdsctools.readthedocs.io/). Microsatellite instability was used as covariate and threshold on minimum size of the positive and negative population for each feature set to 3. The threshold for significance of the *P* value was set to 0.05, but no correction for multiple hypothesis testing was applied.

Kruskal–Wallis rank sum test (one-way ANOVA on ranks) was used to test whether estimated parameters (from bootstrap) derive from the same distribution for all 14 cell lines (null hypothesis rejected if different for at least one group). Effect size *w* was computed as
where
is the statistics from the test and *N* is the number of observations. Effect size >0.5 is considered large. *P* values were Bonferroni corrected to compensate for multiple comparisons, and threshold set to 0.05. Wilcoxon rank sum test was then used for *post hoc* pairwise test on parameters passing Kruskal–Wallis rank sum test, and effect size
, where *Z* is the statistics from the test and *N* is the number of observations. Effect size >0.5 is considered large. *P* values were Bonferroni corrected and threshold was set to 0.05. Rank type tests were preferred over parametric tests because they are highly robust against non-normality.

### Drug combination sensitivity experiments

Experiments were performed using MEKi trametinib (0.005 μmol/L) as anchor drug, testing its effect on our 14 colorectal cancer cell lines alone and in combination with 5 concentrations (from 10 μmol/L with 4-fold dilutions: 0.039 μmol/L, 0.156 μmol/L, 0.625 μmol/L, 2.5 μmol/L, 10 μmol/L) of two GSKi (SB 216763 and CHIR-99021). For available replicates, median value was computed. Dose–response curves were normalized between 0 and 100 with 100 corresponding to the negative control (cells with no treatment) and 0 to the positive control (cells treated with staurosporine 2 μmol/L). Same DMSO volume was maintained for all experiments (negative and positive controls, single drug, and drug combination). Cell number was measured after 3 days of constant drug exposure using CellTiter-Glo reagent as described by the manufacturer (Promega).

## Results

To identify biomarkers of drug response from dynamic logic models, we proceeded as follows (Fig. 1). First, for 14 colorectal cancer cell lines, we systematically obtained antibody-based signaling data on how key molecules respond to perturbations in the signaling network. By combining these data with a literature-derived prior knowledge network, we generated cell line–specific computational models for each of the 14 cell lines. Each model consists of logic ordinary differential equations (22) representing the dynamics of the signaling network, which was parameterized and refined by fitting it to the experimental data using CellNOpt (23), as illustrated in Fig. 1A. Second, we retrieved drug sensitivity data for these 14 cell lines from the GDSC (21), filtering for drugs targeting kinases acting either inside the modeled network or as direct neighbors (Fig. 1B). Third, we investigated correlations between model parameters and drug response to select strong associations using Elastic Net regression (25), as illustrated in Fig. 1C. All steps will be detailed in the following sections.

### Systematic perturbation data for signaling in a large colon cancer cell line panel

To generate a solid dataset to build cell line–specific signaling models and characterize signaling heterogeneity in colon cancer, we performed a large-scale signaling perturbation screen in a panel of 14 genetically heterogeneous colorectal cancer cell lines. We chose the cell lines to reflect the heterogeneity of colon cancer (26), covering typical recurrent mutations in colorectal cancer (such as *KRAS, TP53, APC, BRAF*). Four cell lines had microsatellite instabilities (MSI), and our cell lines covered all classes of the consensus molecular subtype (CMS) described for colorectal cancer patients (27) and computed for colorectal cancer cell lines (see Fig. 2; ref. 28).

Our perturbation data consisted of changes in 14 phosphoproteins measured after the signaling network was systematically perturbed in 43 different scenarios. Perturbations consisted of combinations of 7 small-molecule inhibitors (targeting the kinases IKK, MEK, PI3K, BRAF, TAK1, TBK1, mTOR) and 5 ligand-stimulating receptors (EGF, HGF, IGF1, TGFβ, TNFα; see details in Materials and Methods). Altogether, we obtained 8,428 different experimental data points. When visualizing the perturbation-response profiles (Fig. 2), we noticed that they differ strongly between cell lines. For example, several cell lines showed strong attenuation of AKT signaling upon mTOR and PI3K inhibition (blue areas in the top right of each cell line, such as for CAR1, HCT116, and DIFI), while others did not (such as SW837 and SW620). Notably, these differences were not correlated with consensus subtype or mutation status of, for example, *PIK3CA*. Cluster analyses of the perturbation-response profiles also unveil no clear correlation between genotype data or consensus subtype and perturbation response data, highlighting the idea that the state and responsiveness of the signaling is not fully determined by mutations, and that the obtained dataset provides information that is not contained in genomic data.

We next aimed to build computational network models from the perturbation-response signaling data, to get more insights into the underlying changes in the signaling networks that lead to these heterogeneous responses.

### Compiling a signaling network as a basis for the computational model

As a basis to model the signaling network, we compiled a comprehensive prior knowledge network (PKN; shown in Supplementary Fig. S4) including our measured and perturbed pathways and their connections, using literature and OmniPath (29), a curated ensemble of multiple pathway resources. Overall, the network had 49 nodes (signaling molecules) and 86 edges (regulatory interactions) and could be interpreted as a logic model. As the PKN contained many nodes that were neither measured nor perturbed, we compressed it, as described in ref. 30, to reduce model complexity without affecting logic consistency, while maintaining interactions between measured or perturbed nodes. The model definition considered that PLX4720, a selective BRAF inhibitor, inhibits mutated *BRAF* (V600E), but induces paradoxical activation in wild-type *BRAF* cells (31).

Using the final network, we created a mathematical model that represents the dynamics of the signaling network using logic ordinary differential equations (22). For each node in the network, this model has one parameter *τ _{i}* representing the responsiveness of node

*i*, and for each interaction, there is one parameter

*k*describing the strength of the edge from node

_{j,i}*j*to node

*i*(see Materials and Methods). Note that this model represents a rather generic signaling network, and it is likely that not all interactions are functional in a specific cell line.

### Strategy for model optimization: an *in silico* example

After having established a generic computational model for signaling dynamics, we next aimed to adjust the parameters for each of the 14 cell lines so that we obtain 14 cell line–specific models. A major challenge when attempting to parameterize such large-scale networks is to avoid overfitting, that is, to estimate too many parameters from insufficient data. Building on refs. 23 and 32, we therefore developed a regularization scheme that favors sparser networks (i.e., networks that contain fewer edges and nodes). This is done by adding to the objective function the sum of the absolute parameter values multiplied by a parameter λ (as described in Materials and Methods and Supplementary Material). This way, interactions or responses of nodes that are not supported by the data will be removed and the network will be pruned. The approach is illustrated and characterized in Fig. 3 using an *in silico* model, where the PKN includes 9 nodes and 14 edges. To simulate that some interactions are present in a specific cell line, while others might not, we chose 10 edges (shown in black in Fig. 3A) to be present, and 4 to be absent to generate *in silico* data for 2 readouts (nodes shown in blue in Fig. 3A) recorded after perturbation with 2 stimuli (green nodes) and 2 inhibitors (targeting red nodes). *In silico* data contained in total 10 different perturbation scenarios and were generated for 10 random sets of parameters (one example shown in Fig. 3B).

When we did not include a penalty for complexity (corresponding to λ = 0), the fit of the optimized models to the data was perfect (Fig. 3C with λ = 0); however, parameters could not be well estimated due to model redundancy and low identifiability (Fig. 3D). However, when we increased the penalty for the node parameters *τ _{i}* using the regularization parameter λ (see Materials and Methods), we tended to have worse fit to the data (higher fit error

*Q*) and sparser models (lower sum of estimated parameters

_{LS}*P*), as shown in Fig. 3C. Good values of λ are those on the elbow of the L-shaped curve in Fig. 3C, corresponding to the best compromise between good fit and sparse model. Interestingly, these values (especially λ = 0.001, which represents the most conservative choice in terms of regularization) also correspond to the best accuracy in estimating the parameters in Fig. 3D. Looking at the estimated parameters we could verify that for λ = 0.001 the parameters

*τ*and

_{D}*τ*, corresponding to the unconnected nodes

_{E}*D*and

*E*, were indeed set to zero thanks to the regularization term in the objective function.

### Generation of cell line–specific models

Having established a modeling framework that can handle the complexity of our data and network (Fig. 4A), we generated cell line–specific models for the 14 cell lines, performing bootstrapping also to obtain confidence intervals for the edge parameters *k _{j,i}* (see details in Materials and Methods and Supplementary Methods). The resulting models were comparably sparse across cell lines, with a percentage of values set to zero ranging between 12% and 28% for node parameters

*τ*and 51% and 73% for edge parameters

_{i}*k*. The optimized models showed good fit to the experimental data for all cell lines, with mean squared error (MSE) ranging between 0.013 and 0.036 (see Fig. 4B). Estimated median values for the edge parameters

_{j,i}*k*and node parameter values

_{j,i}*τ*are shown in Fig. 4C and reported in Supplementary Table S2 for all cell lines. Cluster analysis of the model parameters reemphasized that cell lines with similar genetic alterations are not associated with specific states of their signaling networks represented by the parameters (Fig. 4C).

_{i}### Heterogeneity of the cell line–specific models

To further quantify how much the signaling networks differ between cell lines, we statistically investigated the heterogeneity of the estimated edge parameters *k _{j,i}*. Forty-six of the 63 parameters showed a difference in at least one cell line (Kruskal–Wallis rank sum test, adjusted

*P*value <0.05, effect size >0.5, see Materials and Methods). Notably, all other 17 parameters were zero for all cell lines, indicating signaling interactions that are not functional in any of the colon cancer cell lines.

As an attempt to quantify heterogeneity, we tested how often each of the 46 parameters that show variability between cell lines, differs between pairs of cell lines (Wilcoxon rank sum *post hoc* test for all 91 combinations of 14 cell lines, adjusted *P* value < 0.05 and effect size >0.5, see Materials and Methods). The resulting percentages of pairwise differences are visualized as line width and numbers in Fig. 4A, while edge color intensity maps its mean value across cell lines. It is interesting to see that the major signaling routes are among the most variable between cells (PI3K/AKT and BRAF/MEK/ERK); however, also other often overlooked interactions such as feedbacks from RPS6 to PI3K or ERK to BRAF have high variability. The highly variable interaction between IKK and ERK may reflect variability of off-target effects of the chosen IKK inhibitor, as shown previously (8). Importantly for further analysis, we noted that only 34 of 46 edge parameters *k _{j,i}* were different from zero in at least three cell lines, and decided to concentrate our further analysis on these highly heterogeneous parameters.

When investigating the parameters for the nodes, *τ _{i}*, only 1 of 25 parameters was zero across all cell lines and 21 parameters were different from zero in at least three cell lines. Among these, 12 were defined as highly heterogeneous, showing reasonable variability across cell lines (SD > 0.04). From our analysis, we conclude that cell lines are extremely heterogeneous and hence it is important to characterize each cell line with a specific signaling model.

### Association of model parameters with drug response data

We next investigated whether the highly heterogeneous model parameters can be associated with the efficacy of drugs. We reasoned that as many kinase inhibitors target the pathways that show highly variable parameters, these parameters might, in turn, determine the efficiency of these drugs. Drug sensitivity of our 14 colorectal cancer cell lines has been characterized for hundreds of drugs as part of the GDSC project (21). From that dataset, we selected those drugs that target molecules within our network, or molecules that are directly interacting with our network (based on Omnipath; ref. 29). We excluded drugs with low variability in the response across our 14 cell lines (less than three sensitive and three resistant cell lines based on the classification used in ref. 21) and those for which sensitivity data were missing for more than three of our 14 cell lines. For the remaining 27 drugs, we investigated whether the sensitivity of our cell lines (as measured by the IC_{50}) is correlated with any of the model parameters (using cross-validated Elastic Net, details in Materials and Methods). Out of all 1,242 potential correlations, we found in total 146 significant and robust parameter–drug associations. These associations covered 14 of our 27 drugs (Supplementary Figs. S5 and S6). We then performed a similar analysis to obtain associations between drugs and genomic alterations (functional mutations and copy number alterations) acting on PKN nodes or first neighbor (using ANOVA as described in ref. 21, without correction for multiple hypothesis testing due to too large number of comparisons), we obtained associations for 12 of our 27 drugs, five of which also showed associations with model parameters (significant associations are reported in Supplementary Fig. S7).

### Interpretation and validation of the dynamic biomarkers

We then focused on the 9 drugs whose efficacy was not significantly associated with mutations or copy number alterations but instead with pathway biomarkers, to explore how our approach provides insights in otherwise unexplained drug efficacy. Top associations are illustrated in Fig. 5, visualized on the network (PKN) to show, for each drug, both the target(s) (indicated with pills) and the associated pathway parameter(s) (colored edges).

For example, multiple parameters related with MEK–ERK pathway showed association with EGFR signaling inhibitor (HG-5-88-01) and with BRAF inhibitors (SB590885 and dabrafenib). Interestingly, many recent studies reported synergistic effects when combining MEK inhibitors with EGFR inhibitors (8, 33, 34) or with BRAF inhibitors (8, 35, 36) in different cancer types. In particular, focusing on colorectal cancer, combination of EGFR and MEK inhibitors was suggested to overcome resistance to MEK inhibitors in BRAF and KRAS mutants (8, 37) and it was suggested in ref. 34 as a potential therapy for colorectal cancer patients who become refractory to anti-EGFR therapies. A significant increase in response was also shown on a phase I study for colorectal cancer patients treated with a combination of three drugs targeting MEK, BRAF, and EGFR respectively, in BRAF V600E mutants (38) and is currently in phase I/II study (NCT01750918 in clinicaltrials.gov). These examples underline the potential of our approach in suggesting targeted combinatorial therapies.

Another interesting example is the association of different MEK inhibitors (AZD6244, RDEA119, trametinib, PD-0325901, CI-1040) with multiple model parameters related to GSK3 (*τ _{GSK3}, k_{RSKp90,GSK3}, k_{AKT,GSK3}*), as shown in Fig. 6A–C, with each parameter associated with at least two of the inhibitors. The top association among these, between AZD6244 and

*τ*, is shown in the center of Fig. 6D. Our interpretation of these associations is that inhibition of GSK3 would improve the sensitivity for cell lines resistant to MEK inhibitors (i.e., high IC

_{GSK3}_{50}) and with functional GSK3 (i.e., high

*τ*). On the contrary, cell lines with nonfunctional GSK3 (i.e., low

_{GSK3}*τ*) should not be affected by GSK3 inhibitors. To experimentally test this hypothesis, we combined a MEK inhibitor (Trametinib) with two GSK3 inhibitors (SB216763 and CHIR-99021) at increasing concentrations (Materials and Methods) as shown in the external plots in Fig. 6D. As expected, cell lines with nonfunctional GSK3 do not show an increased sensitivity when treated in combination with any of the two tested GSK3 inhibitors, regardless of their response to MEK inhibitors (e.g., HT29 and SNUC2B for a sensitive and a more resistant case, respectively). On the contrary, cell lines with more functional GSK3 tend to show an improved sensitivity under the combinatorial treatment (see HT115, SW1116, SW1463, and CCK81) even at the low drug concentrations tested.

_{GSK3}## Discussion

It is broadly accepted that alterations in signaling pathways largely determine the efficacy of kinase inhibitors used in the clinic. Yet, in only few examples we have a clear understanding of how changes in signaling networks shape drug response. In this study, we systematically investigated how the dynamics of signaling pathways can determine the efficacy of targeted drug treatments and found several examples where quantitative differences in signaling networks are associated with drug sensitivity that cannot be explained by genomic data. In addition, we show that this information can be used to guide combinatorial therapies.

To uncover cell-specific changes in signaling, we constructed mechanistic models of signal transduction for 14 colorectal cancer cell lines based on the changes in phosphorylation in key signaling molecules upon perturbations. We built our models based on differential equations to capture the continuous aspects of signal strength and time dynamics, but using a logic formalism that kept the model simple and with easily interpretable parameters. We could then study how these model parameters, which determine the dynamic behavior of the pathway, are related with the cellular sensitivity to cancer drugs. We found 146 strong associations between model parameters and drug sensitivity, including drugs without genetic markers. This suggests that differential wiring of these pathways is related to the efficacy of the drugs, and these differences could hint at alternative markers.

A key aspect of our study is that we capture and catalog the heterogeneity of signaling dynamics (which is described by the model parameters) for a large panel of colorectal cancer cell lines. It uses a unique dataset, consisting of 14 phospho-proteins measured under 43 differently perturbed conditions, which allowed us to generate and characterize cell line–specific models for 14 cell lines. In our models trained on this data, we found that about half of the model parameters are highly heterogeneous in our panel of cell lines and that this heterogeneity cannot be explained by genomic alterations alone. As these findings were obtained with a rather focused set of antibody-based assays, we expect that many more associations can be found by a broader characterization. Mass spectrometry–based approaches, while still more costly and laborious than antibody-based methods, are under active development and likely to enable such an analysis in the near future (39).

Associations between model parameters and drug sensitivity were used to define interesting points of intervention to overcome resistance to specific drugs in specific cell lines, thus paving the way for personalized combinatorial therapies. Some of the combinatorial therapies suggested using this approach, like the combination of MEK inhibitors with EGFR and/or BRAF inhibitors, are supported by literature (as reported in the Results) and are currently in clinical trials. We found a rather surprising association between GSK3-related parameters and the efficacy of MEK inhibitors, which suggests that GSK3 inhibition would sensitize cells for drugs targeting MEK, which we validated experimentally. This interaction could not be unveiled using genomic data, but required information about the dynamics of signaling. Interestingly, the inhibition of GSK3 showed promising anticancer effect in recent studies on cancer cell lines (40, 41) and its potential interaction with the MEK–ERK pathway has also been suggested in this context (41, 42) but, as far as we know, no proven synergistic effect had previously been reported.

Our study clearly highlights the influence of intertumor heterogeneity in signaling networks on drug response; however, each colorectal tumor itself is also a very heterogeneous system. The tumor cells form hierarchies with stem cells and more differentiated cells, giving rise to a mix of different interacting cell types, including fibroblasts and immune cells (43, 44). As each of these cell types has different signaling capabilities they will likely respond differently to the targeted drugs, and in turn the drugs will influence the heterocellular composition of the tumors. Therefore, it will be very important to study signaling response in these heterogeneous cells in future, such as by single-cell technologies, to fully appreciate the complexity of the system *in vivo*.

Our findings underscore the value of studying the dynamics of signaling pathways to better understand tumor phenotypes and to exploit this knowledge to suggest new therapeutic strategies (45). Such an approach is fundamentally different from the currently common characterization of various “omics” layers at the basal level and it can be exploited to anticipate the dynamic adaptation mechanisms underlying drug resistance (17). Accordingly, we were able to find various pathway dynamic biomarkers for drugs with no genetic biomarkers. Hence, we believe that a perturbation-based strategy, even if restricted to fewer genetic backgrounds and only monitoring selected proteins, can provide a complementary strategy for biomarker discovery in cancer and beyond.

## Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

## Authors' Contributions

**Conception and design:** F. Eduati, B. Klinger, N. Blüthgen, J. Saez-Rodriguez

**Development of methodology:** F. Eduati, V. Doldàn-Martelli, T. Cokelaer, J. Saez-Rodriguez

**Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.):** A. Sieber, M.J. Garnett

**Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis):** F. Eduati, V. Doldàn-Martelli, B. Klinger, T. Cokelaer, M. Dorel, N. Blüthgen, J. Saez-Rodriguez

**Writing, review, and/or revision of the manuscript:** F. Eduati, B. Klinger, M.J. Garnett, N. Blüthgen, J. Saez-Rodriguez

**Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases):** T. Cokelaer, F. Kogera, M. Dorel

**Study supervision:** N. Blüthgen, J. Saez-Rodriguez

## Grant Support

N. Blüthgen was supported by the Berlin Institute of Health and by the Federal Ministry of Research and Education (BMBF; grants OncoPath and ColoSys). F. Eduati was supported by European Molecular Biology Laboratory Interdisciplinary Post-Docs (EMBL EIPOD) and Marie Curie Actions (COFUND).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked *advertisement* in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

## Acknowledgments

We thank Markus Morkel and members of the Saez laboratory for useful feedback, especially Luz Garcia-Alonso for help on mutations and Attila Gabor for fruitful discussion on model definition and David Henriques for support with CNORode. We thank Howard Lightfoot for help and assistance with drug combination data.

## Footnotes

**Note:**Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

- Received January 10, 2017.
- Revision received March 17, 2017.
- Accepted March 31, 2017.

- ©2017 American Association for Cancer Research.