Abstract
Resistance to targeted cancer therapies is an important clinical problem. The discovery of antiresistance drug combinations is challenging as resistance can arise by diverse escape mechanisms. To address this challenge, we improved and applied the experimentalcomputational perturbation biology method. Using statistical inference, we build network models from highthroughput measurements of molecular and phenotypic responses to combinatorial targeted perturbations. The models are computationally executed to predict the effects of thousands of untested perturbations. In RAFinhibitor resistant melanoma cells, we measured 143 proteomic/phenotypic entities under 89 perturbation conditions and predicted cMyc as an effective therapeutic cotarget with BRAF or MEK. Experiments using the BET bromodomain inhibitor JQ1 affecting the level of cMyc protein and protein kinase inhibitors targeting the ERK pathway confirmed the prediction. In conclusion, we propose an anticancer strategy of cotargeting a specific upstream alteration and a general downstream point of vulnerability to prevent or overcome resistance to targeted drugs.
https://doi.org/10.7554/eLife.04640.001eLife digest
Drugs that target the activity of specific genes could potentially form precise cancer therapies. Some cancers, including the aggressive skin cancer called melanoma, initially respond well to such treatments. However, resistance to drugs develops quickly, leading to the rapid regrowth of the tumors.
Resistance can develop in a number of ways. For example, to prevent the drug from working or to compensate for the effects of a drug, cancer cells can adapt their signaling processes or acquire genetic mutations or other modifications that affect how genes are expressed.
A welldesigned combination of drugs that targets multiple molecular pathways can make it harder for cells to resist treatment, as this limits the number of available ‘escape’ pathways that bypass the drug targets. However, it is difficult to accurately predict how a cell will respond when treated with a particular drug, making it extremely challenging to design effective drug combinations.
In 2013, researchers developed a technique to build predictive models of cellular response pathways based on data collected from perturbation experiments followed by mathematical modeling. Now, Korkut et al.—including several of the researchers involved in the 2013 work—have refined this technology and applied it to the problem of preventing drug resistance in cancer cells.
Computer simulations that used the mathematical models suggested a particular strategy of ‘upstream–downstream targeting’ in cells that were insensitive to the clinically successful drug vemurafenib (an inhibitor of RAF proteins, which are often mutated in cancers). In the landscape of signaling pathways, the target of the upstream drug is on or near the mutated RAF protein. cMyc, the indirect target of the downstream drug helps to express genes that trigger signals that cause the cells to grow. Inhibiting both targets in parallel may have the dual advantage of blocking the activation of the tumorspecific growth pathway while reducing the cancer cells' attempts to bypass the activation block.
An initial test of the designed drug combination required moving from computer simulations to the laboratory using cell cultures originally derived from melanoma tumors. When Korkut et al. applied the drug combination, the combined treatment successfully blocked cell growth. The results suggest that the datadriven computer modeling strategy termed perturbation biology could be a useful tool for identifying effective cancer drug combinations for further preclinical research, possibly followed by clinical trials.
https://doi.org/10.7554/eLife.04640.002Introduction
Drug resistance in cancer treatment
The inhibition of key oncogenes with targetspecific agents elicits dramatic initial response in some cancers such as prostate cancer and melanoma (Bollag et al., 2010; Clegg et al., 2012). Even for the most successful singleagent targeted therapies, however, drug resistance eventually emerges leading to rapid progression of metastatic disease (Garraway and Janne, 2012). The mechanism of drug resistance may involve selection of resistant subclones, emergence of additional genomic alterations, and compensating interactions between alternative signaling pathways (Choi et al., 2007; Huang et al., 2011; Wagle et al., 2013). One potential solution to overcome drug resistance is to combine targeted drugs to block potential escape routes (Fitzgerald et al., 2006). Therefore, there is currently a need for systematic strategies to develop effective drug combinations.
Drug resistance in melanoma
Targeted therapy has been particularly successful in treatment of melanoma. BRAFV600E gainoffunction mutation is observed in ∼50% of melanomas (Davies et al., 2002). Direct inhibition of BRAFV600E by the RAF inhibitor (RAFi) vemurafenib and inhibition of the downstream MEK and ERK kinases have yielded response rates of more than 50% in melanoma patients with this mutation (Chapman et al., 2011; Flaherty et al., 2012b). At the cellular level, inhibition of the ERK pathway leads to changes in expression of a set of critical cell cycle genes (e.g., CCND1, MYC, FOS) and feedback inhibitors of ERK signaling (e.g., DUSP, SPRY2) (Pratilas et al., 2009). Resistance to vemurafenib emerges after a period of ∼7 months in tumors that initially responded to singleagent therapy (Chapman et al., 2011; Sosman et al., 2012). Multiple RAFi and MEKi (e.g., PD0325901, Trametinib) resistance mechanisms, which may involve alterations in NRAS/ERK pathway (e.g., NRAS mutations, switching between RAF isoforms) or parallel pathways (e.g., PTEN loss), have been discovered in melanoma (Johannessen et al., 2010; Nazarian et al., 2010; Poulikakos et al., 2010; Xing et al., 2012).
The alterations associated with drug resistance may preexist alone, in combinations, or emerge sequentially and vary substantially between patients (Van Allen et al., 2014). Effective drug combinations may target diverse resistance mechanisms. Despite anecdotal success, conventional methods and combinatorial drug screens generally fail to come up with effective combinations due to the genomic complexity and heterogeneity of tumors (Zhao et al., 2014). In order to more effectively nominate drug combinations, we propose to employ systemwide models that cover interactions between tens to hundreds of signaling entities and can describe and predict cellular response to multiple interventions. There have been prior attempts to construct such signaling models. De novo and datadriven quantitative models were able to cover only a few signaling interactions and therefore had limited predictive power (Nelander et al., 2008; Bender et al., 2011; Klinger et al., 2013; Oates et al., 2014). Qualitative or discrete models can cover more interactions but typically lack the capability of generating quantitative predictions (SaezRodriguez et al., 2009; Breitkreutz et al., 2010; SaezRodriguez et al., 2011). Detailed physicochemical models derived using generic biochemical kinetics data can be quite comprehensive and quantitative but typically lack sufficient celltype specificity required for translationally useful predictions (Chen et al., 2009).
We construct comprehensive, celltype specific signaling models that quantitatively link drug perturbations, (phospho)proteomic changes, and phenotypic outcomes (Figure 1). The models capture diverse signaling events and predict cellular response to previously untested combinatorial interventions. In order to generate the training data for network modeling, we first perform systematic perturbation experiments in cancer cells with targeted agents. Next, we profile proteomic and phenotypic response of cells to the perturbations. The celltype specific response data serve as the input for network inference. In this study, we also incorporate prior pathway information from signaling databases to narrow the parameter search space and improve the accuracy of the models. For this purpose, we developed a computational tool (Pathway Extraction and Reduction Algorithm [PERA]) that automatically extracts priors from the Pathway Commons signaling information resource (Demir et al., 2010; Cerami et al., 2011).
Even in the presence of large training data and priors, network inference is a difficult problem due to the combinatorial complexity (i.e., exponential expansion of the parameter search space with linear increase of parameters). For example, to infer a network model with 100 nodes using Monte Carlobased methods, we in principle would need to cover a search space that includes ∼2^{(100 × 100)} network models—a computationally impossible task. To circumvent this problem, we previously had developed a network modeling algorithm based on belief propagation (BP), which replaces exhaustive onebyone searches over many individual network models by a search over probability distributions representing sets of lowerror network models (Molinelli, Korkut, Wang et al., 2013). The algorithm enables us to construct models that can predict response of hundreds of signaling entities to any perturbations in the space of modeled components. Here, we improved the perturbation biology method through automated incorporation of prior information (signaling interactions from databases) to obtain more accurate network models. The prior information is extracted from the signaling databases using the newly developed PERA tool. To improve predictive power and preserve celltype specificity, we use prior information as soft restraints on search space through use of a probabilistic error model of priors, that is, the algorithm rejects interactions that do not conform to the experimental training data and predicts novel interactions not sampled in the priors (see Equations 2–7). To derive richer and more informative network models, we also scaled up the implementation to deal with a larger number of proteomic and phenotypic nodes (see ‘Materials and methods’ for details and comments on potential future opportunities in perturbation biology). To quantitatively predict cellular response to combinatorial perturbations, we simulate the fully parameterized network models with in silico perturbations until the system reaches steady state (Figure 1). The steadystate readout for each proteomic and phenotypic entity (i.e., system variables) is the predicted response to the perturbations.
In this study, we improved and applied the perturbation biology method to devise a potentially generalizable strategy for overcoming resistance to targeted cancer therapies. We constructed celltype specific network models of signaling from perturbation experiments in RAFiresistant melanoma cells (SkMel133 cell line). The melanoma cells used for network modeling have a BRAFV600E mutation and homozygous deletions in PTEN and CDKN2A. The models quantitatively link 82 (phospho)proteomic nodes (i.e., molecular concentrations) and 12 protein activity nodes with 5 cellular phenotype nodes (e.g., cell viability). As shown by crossvalidation calculations, use of prior information significantly improved the predictive power of the models. Once the predictive power was established, we expanded the extent of the drug response information from a few thousand experimental data points to millions of predicted node values. Based on the predictions, which cannot be trivially deduced from experimental data, we nominated cotargeting of cMyc and BRAF or MEK as a potential strategy to overcome RAFiresistance. To test the prediction, we first experimentally showed that the BET bromodomain inhibitor, JQ1, reduces cMyc expression. Next, we showed that JQ1 blocks the growth of RAFiresistant SkMel133 cells in synergy with RAF/MEK signaling inhibition, and in this context, overcomes the drug resistance. Based on these results, we put forward the falsifiable hypothesis that cotargeting a specific upstream alteration and a general downstream point of vulnerability is a good strategy to prevent or overcome resistance to targeted drugs.
Results
Experimental response maps from proteomic/phenotypic profiling
Drug perturbation experiments in melanoma cells
We performed systematic perturbation experiments in malignant melanoma cells (Figure 2A) to generate a rich training set for network inference. The RAFiresistant melanoma cell line SkMel133 was treated with combinations of 12 targeted drugs (Figure 2A, see Supplementary file 1A for drug targets and doses, Supplementary file 1C for all perturbation conditions). The perturbations consisted of systematic paired combinations and multiple doses of single agents. This procedure generated 89 unique perturbation conditions, which targeted specific pathways including those important for melanoma tumorigenesis such as ERK and PI3K/AKT (Haluska et al., 2006).
Proteomic/phenotypic profiles
A key aspect of the data acquisition for network inference is combining the proteomic and cellular phenotypic data so that the resulting models quantitatively link the proteomic changes to global cellular responses. Toward this objective, we profiled the melanoma cells for their proteomic and phenotypic response under 89 perturbation conditions (Figure 2B,C). We used reverse phase protein arrays (RPPA) to collect drug response data for 138 proteomic (total and phosphoprotein levels) entities in all conditions (Tibes et al., 2006). In parallel, we measured phenotypic responses, including cell viability and cell cycle progression (i.e., G1, S, G2, G2M arrest phenotypes) in all conditions (Figure 2B).
The response map
The highthroughput phenotypic and proteomic profiles form a response map of cells to systematic perturbations (Figure 2). The response map provides contextspecific experimental information on the associations between multiple system variables (i.e., proteomic entities) and outputs (i.e., phenotypes) under multiple conditions (i.e., perturbations). We demonstrate through hierarchical clustering of the map that each targeted drug induces a distinct proteomic response, and drugs targeting the same pathway lead to overlapping responses in the SkMel133 cells (Figure 2B). Through a clusteringdriven pathway analysis, we further show that functionally related proteins (i.e., proteins on same or related pathways) respond similarly to targeted agents (Figure 2C, Figure 2—figure supplements 1–3).
The response map can be described as a set of uncoupled pairwise associations between proteomic and phenotypic entities. We provide an example of an association network model (Figure 2—figure supplement 4). The association network was inferred using the experimental response map (Figure 2) and graphical Gaussian models, which utilizes partial correlations for network inference (Schafer and Strimmer, 2005). Such models, however, are not sufficient for achieving systematic predictions as they do not capture the nonlinear nature of the couplings between the entities and cannot be executed with in silico perturbations. Therefore, we built quantitative models using the experimental response map. The resulting models describe the coupled nature of the interactions between proteins and cellular events, as well as the nonlinear dynamics of cellular responses to drug perturbations.
Quantitative and predictive network models of signaling
Network models
Next, we used the experimental response map (Figure 2) and the BPbased inference strategy (Figure 1) to build quantitative network models of signaling in melanoma. In the models, each node quantifies the relative response of a proteomic or phenotypic entity to perturbations with respect to the basal condition. Consequently, proteomic entities that do not respond to even a single perturbation condition, do not contribute any constraints for inference. We eliminated such entities from the network modeling with a signaltonoise analysis and included 82 of the 138 proteomic measurements in the modeling (see ‘Materials and methods’, Supplementary file 1B). In addition to the proteomic nodes, the models contained 5 phenotypic nodes and 12 ‘activity nodes’, which represent the 12 drugs and couple the effects of the targeted perturbations to the other nodes in the network. In total, network models contained 99 nodes. BP algorithm generates the probability distribution of edge strengths for every possible interaction between the nodes. The BPguided decimation algorithm (see ‘Materials and methods’, Figure 1—figure supplement 1) instantiates distinct network model configurations from the probability model (Montanari et al., 2007).
The mathematical formulation of the BPbased network inference is suitable for both de novo modeling (i.e., modeling with no prior information) and modeling using prior information (see ‘Materials and methods’). Here, we used prior information to infer models with higher accuracy and predictive power compared to de novo models. Using the PERA computational tool, we derived a generic prior information model from Reactome and NCINature PID databases, which were stored in Pathway Commons (Cerami et al., 2011). The prior information network contains 154 interactions spanning multiple pathways (Figure 4—figure supplement 1). Next, we added a prior prize term to the error model (Equation 2) to restrain the search space by favoring the interactions in the prior model. It is critical that the prior information does not overly restrain the inferred models and the algorithm can reject incorrect priors. To address this problem, we inferred network models using the pathway driven and randomly generated prior restraints. The statistical comparison of the networks inferred with actual (i.e., reported in databases) and random prior models indicated that the inference algorithm rejected a significantly higher number of prior interactions when randomly generated priors were used for modeling (Figure 4—figure supplement 3). Finally, we integrated the experimental data and prior information to generate 4000 distinct and executable model solutions with low errors (i.e., sets of model parameters, W_{ij}, for which the model equations best reproduce the experimental response map) using the BPbased strategy.
Use of prior information improves the predictive power of models
Cross validations with and without prior information
We addressed the question whether BPderived models have predictive power and whether use of prior information introduces further improvement. To assess the predictive power of the network models (i.e., predicting the response to untested perturbations), we performed a leavekout cross validation (Figure 3A). In each validation calculation, we withheld the response profile to every combination of a particular drug (e.g., RAFi) and all other drugs (leave11out cross validation) (Figure 3B, Figure 3—figure supplement 1). For each crossvalidation calculation, this procedure created a partial training data set that contains response to combinations of 11 drugs and two different doses of a single drug totaling to 78 unique conditions (Supplementary file 1B,C).
First, we constructed both de novo (i.e., without any prior information) and prior information guided network models using two partial data sets from which responses to combinations of drugs with RAFi or AKTi were withheld. Next, we predicted the response by executing the models with in silico perturbations that corresponded to the withheld experimental conditions. Finally, we compared the hidden and predicted response data from models generated de novo or with prior information.
Restraining inference with prior information improves the predictive power of models
The comparison between the predicted and the withheld experimental profiles suggests that the de novo network models have considerable predictive power, and the use of prior information in modeling introduces significant improvement in the prediction quality (Figure 3B). Use of prior information increased the cumulative correlation coefficient between predicted and experimental response data from 0.72 to 0.84 and from 0.71 to 0.81 for RAFi and AKTi, respectively. The demonstrated predictive power of the models suggests that the models are suitable for systematically predicting response to perturbation combinations not sampled in the training set and generating testable hypotheses that link external perturbations (e.g., targeted drugs) to cellular response.
Network models identify contextdependent oncogenic signaling in melanoma
Network modeling and the average model
We generated quantitative network models with the complete experimental response profile and prior information to investigate oncogenic signaling in melanoma. The resulting network models resemble conventional pathway representations facilitating their comparison with the biological literature (Le Novere et al., 2009), but the interaction edges do not necessarily represent physical interactions between connected nodes. Analysis of the ensemble of network model solutions reveals that a set of strong interactions is shared by a majority of the inferred network models. On the other hand, some interactions have nonzero edge strength (W_{ij} ≠ 0) values only in a fraction of the models (see Figure 4—figure supplement 2 for analysis of the edge distribution in models). As a first step of detailed analysis and for the purpose of intuitive interpretation, we computed an average network model (Figure 4A), which is obtained by averaging the interaction strength (W_{ij}) for each node pair, ij, over all individual model solutions. This average network model highlights the interactions with high W_{ij} values if shared by the majority of the distinct solutions. Although the average model cannot be simulated to predict system response, it is particularly useful for qualitative analysis of the inferred signaling interactions.
Global analysis of average models
The average network model provides a detailed overview of the signaling events in melanoma cells (Figure 3A). The average model contains 202 unique interactions (127 activating and 75 inhibitory interactions) between 99 signaling entities. 89 of the 154 interactions in the prior model conform to the experimental data, and therefore, are accepted in the majority of the model solutions by the inference algorithm and included in the average model (Figure 4—figure supplement 4). Given that the average model covers interactions from multiple signaling pathways and is more complex than the pathway diagrams presented in most review papers, even the qualitative analysis of the model is highly challenging.
Network models capture known signaling pathways
In order to simplify the analysis of the average model solution, we fragmented the global network diagram into subnetworks (Figure 4). Each subnetwork is a simplified representation of the signaling events in canonical pathways such as those that fall into ERK, PI3K/AKT, and cell cycle pathways (Figure 4C–E). The subnetwork diagrams indicate that models recapitulate many known interactions in pathways, which are important in melanoma tumorigenesis (e.g., PI3K/AKT and ERK) and nominate previously unidentified interactions (see Figure 4 legend). It is, however, not possible to predict the cellular response to untested drug perturbations through qualitative analysis of the inferred interactions. We use quantitative simulations with in silico perturbations to both decode the signaling mechanisms and more importantly systematically predict cellular response to combinatorial drug perturbations.
Combinatorial in silico perturbations generate an expanded proteomic/phenotypic response map
Model execution with in silico perturbations
Thanks to their ordinary differential equation (ODE)based mathematical descriptions, the models can be executed to predict cellular responses to novel perturbations (Nelander et al., 2008). The systematic predictions from the models go beyond the analysis of few particular edges in the system and capture the collective signaling mechanisms of response to drugs. We execute the parameterized model ODEs (Equation 1) with in silico perturbations acting on node (i) as a real numbered u(i) value until all the system variables (i.e., node values, {x_{i}}) reach to steady state (Figure 5A,B).
Prediction of phenotypic responses
The simulations expand the size of the response map by three orders of magnitude (i.e., from few thousand experimental response data to millions of predicted responses) (Figure 5C, Figure 5—figure supplement 1). Once we had the predicted response profiles, we searched for specific perturbations that may induce desired phenotypic changes even when cells are treated with drugs at physiologically relevant doses. For this purpose, we identified the top ranked perturbation conditions in terms of their predicted influence on each phenotype and nominated them for further experimental testing (Table 1, Supplementary file 1F). Not surprisingly, the execution of models quantitatively recaptured the experimentally observed associations between the drug perturbations and the phenotypic responses. For example, targeting PKC or CDK4 with specific kinase inhibitors lead to reduction of cell viability according to the simulations, which can also be directly observed from the experiments. However, CDK4 and PKC inhibitors substantially reduced SkMel133 cell viability only at high doses, as in the original perturbation experiments (>2 μM) (Figure 6—figure supplement 1). More importantly, the simulations allowed us to identify effective perturbation combinations that cannot be trivially deduced from the experimental data (Table 1, Figure 5D–H, Supplementary file 1F).
The execution of the models identified a set of perturbation conditions that are likely to induce a strong phenotypic response. The listed predicted targets (Table 1, Supplementary file 1F) will serve as a useful resource to nominate intervention strategies to overcome RAFi resistance in melanoma. In particular, in silico perturbations of cMyc lead to the strongest G1arrest response when cMyc is cotargeted with BRAF, MEK, and cyclin D1 (Figures 5D–H, 6A). The G1arrest phenotype is of particular interest given that induction of G1arrest is a promising strategy for controlling melanoma growth (Solit et al., 2006) and various genomic alterations associated with G1/S transition are involved in tumor progression (Rother and Jones, 2009). Consistent with the vast literature on cMyc's role in genesis of many cancers (Dang, 2012), predictions indicated that cMyc linked multiple pathways such as ERK and PI3K/AKT to regulation of cell cycle arrest (Table 1; also see Figure 6—figure supplement 2 for the highly complex cMyc response to targeted perturbations in the underlying RPPA data). As neither cMyc nor its direct regulators were inhibited in the perturbation experiments, the predictions from the models were nontrivial and we decided to test these predictions experimentally.
Cotargeting cMyc with MEK or RAF is synergistic in melanoma cells
We predicted through quantitative simulations that melanoma cells were arrested in G1phase of the cell cycle when cMyc was targeted alone or in combination with other proteins, particularly BRAF, MEK, and cyclin D1 (Figure 6A). We experimentally tested the key prediction from the network models. In order to target cMyc expression, we treated melanoma cells with the BET bromodomain inhibitor, JQ1, as a single agent and in combination with MEKi (PD0325901) or RAFi (vemurafenib). JQ1 directly targets bromodomains, especially those of the bromodomain protein BRD4, which marks select genes including MYC on mitotic chromatin. Inhibition of the BRD4 bromodomains with JQ1 downregulates MYC mRNA transcription and leads to G1 cell cycle arrest in diverse tumor types, such as multiple myeloma (Delmore et al., 2011; Loven et al., 2013; Puissant et al., 2013).
First, we asked whether we could affect cMyc levels in SkMel133 cells using JQ1. As measured by Western blot experiments, cMyc protein expression is reduced in response to JQ1 alone. cMyc protein levels are further reduced when the cells are treated with a combination of JQ1 and MEKi or RAFi (Figure 6B).
To directly test the key prediction from the perturbation biology models, we measured the cell cycle progression response of melanoma cells to JQ1 in combination with the RAF and MEK inhibitors. We observed a strong synergistic interaction between JQ1 and RAFi (Figure 6C,D). 51% and 46% of melanoma cells are in G1stage 24 hr after treatment with JQ1 (500 nM) and RAFi (200 nM), respectively, while 39% of cells are in G1stage in the absence of any drug. On the other hand, when cells are treated with the combination of JQ1 and RAFi, a drastic increase in the fraction of cells arrested in G1stage (84%) is observed. The single agent MEKi (50 nM) induces a strong G1arrest phenotype in SkMel133 cells (88% G1stage in MEKitreated cells vs 39% in nondrug treated cells). The combination of MEKi with JQ1 arrests an even higher fraction of the cells (92%) in the G1stage (Figure 6—figure supplement 3).
Before assessing the effect of JQ1MEKi/RAFi combination on viability of melanoma cells (SkMel133), we tested the effect of single agent JQ1 and found that the melanoma cells were considerably sensitive to single agent JQ1 treatment (cell viability IC50 = 200 nM). The sensitivity of SkMel133 to JQ1 is similar to those of A375 and SkMel5 lines (RAFi/MEKi sensitive, carrying BRAFV600E mutation) to another BRD4 inhibitor, MS417 (Segura et al., 2013). The observed sensitivity is also comparable to those of multiple myeloma and MYCNamplified neuroblastoma cell lines, reported to be potentially JQ1sensitive tumor types (Delmore et al., 2011; Puissant et al., 2013), and substantially higher than those of lung adenocarcinoma and MYCNWT neuroblastoma cell lines (Lockwood et al., 2012; Puissant et al., 2013).
We tested the effect of combined targeting of cMyc with MEK or BRAF on cell viability in SkMel133 cells (Figure 6E). Strikingly, when combined with JQ1 (120 nM), cell viability is reduced by 50% with 120 nM of RAFi (PLX4032), whereas the IC50 for single agent RAFi is >1 μM in RAFiresistant SkMel133 cells. Similarly, when combined with 5 nM MEKi (PD901), viability of SkMel133 cells is reduced by 50% with 100 nM of JQ1, an IC50 value, which is close to those of the most sensitive multiple myeloma cell lines (Delmore et al., 2011). At higher doses (IC80), JQ1 is synergistic with both MEKi (combination index, CI_{85} = 0.46) and RAFi (CI_{85} = 0.47) in SkMel133 cells. At intermediate doses, JQ1 synergizes with RAFi (CI_{50} = 0.65) and has near additive interaction with the MEKi (CI_{50} = 0.85) (Figure 6F). Consistent with the observed synergy at high doses, both JQ1 combinations significantly improve the maximal effect level (A_{max}, response to the drugs at highest doses), leading to lower cell viability beyond the levels reached by treatment with any of the agents alone. The observed improvement in A_{max} is particularly important since a subpopulation of cancer cells usually resist treatment even at highest possible drug doses. Treatments with drug combinations, such as those tested here will overcome or delay emergence of drug resistance if they can shrink the size of this resistant subpopulation (i.e., lead to improved A_{max}).
Discussion
We performed a series of systematic perturbation experiments and, measured proteomic and phenotypic responses in each perturbation condition. Using the highthroughput response profiles, we generated network models of signaling in melanoma cells to systematically predict cellular responses to untested drug perturbations. Our modeling algorithm integrates highthroughput drug response profiles and pathway information from signaling databases. The scale and the predictive power of the models were beyond the reach of the previously available network modeling methods, and we made substantial progress in the design of experiments, acquisition of large profiling data sets, and improving a very efficient probabilistic network inference method. Based on the predictions from the datadriven network models, we found that cotargeting MEK or BRAF with cMyc leads to synergistic responses that render a melanoma cell line more sensitive to inhibition of BRAF. Beyond nomination of effective drug combinations in cell lines, we suggest that the perturbation biology method more generally paves the way for modeldriven quantitative cell biology with diverse applications in many fields of biology.
The problem of insensitivity or resistance to targeted therapy is both important and complicated. Oncogenic alterations that affect sensitivity to targeted therapies may preexist in the tumor before treatment or emerge sequentially as the result of treatment. Tumors can escape drug treatment through diverse routes and dominant mechanisms of resistance are the exception (such as the T790M mutation in EGFR, which confers resistance to EGFR inhibitors in many patients as well as in the laboratory [Pao et al., 2005]), so the identification of effective countermeasures is in general difficult. Typically researchers have identified mechanisms of resistance in individual cases and proposed countermeasures based on the particular mechanism. For example, PTEN loss was associated with resistance to RAFi (in BRAFV600E mutated melanoma) and as a result joint inhibition of BRAF and PI3K was proposed as a combination therapy (Villanueva et al., 2010; Xing et al., 2012). Based on similar logic, other plausible therapeutic combinations that target particular escape routes in RAFiresistant melanomas have been put forward, including combinations that target downstream elements of the ERK pathway (e.g., RAFi + MEKi) (Lito, Pratilas et al., 2012; Flaherty et al., 2012a), or upstream receptor tyrosine kinases (Prahallad et al., 2012; Sun et al., 2014). For almost all cases, however, tumor cells rapidly adapt and escape through alternate routes after a brief period of response to the combination. Even the FDAapproved RAFi + MEKi combination has only achieved limited improvements in the median progressionfree survival in BRAFV600E mutated melanoma patients compared to single agent RAFi treatment (Long et al., 2014). A general solution for resistance to the RAFi in BRAFmutated melanoma, in particular, is not yet in sight.
Aiming to develop a more general strategy for preventing the emergence of drug resistance, we have taken a more systematic approach to the discovery of effective combination therapy and cast a wider net by perturbing cancer cells with an entire set of targeted inhibitors and combinations of inhibitors and observing the proteomic signaling response using hundreds of antibodies covering diverse pathways. We expected to generate plausible predictions for a set of effective combinations and to be able to rank their relative importance, to guide validation experiments. The perturbation biology approach used here identified a highly ranked cotarget of RAF, downstream of the ERK pathway, the oncogene cMyc. To test our predictions, we targeted cMyc using the epigenetic drug JQ1, a BRD4 inhibitor that negatively affects cMyc transcription (Delmore et al., 2011). We treated SkMel133 cells with combinations of JQ1 and MEKi or RAFi and showed that both combinations synergistically induced G1arrest and decreased cell viability.
A key criterion of pharmacological efficacy is the ability of a drug to induce cellular responses at doses sufficient to inhibit its immediate molecular targets. In SkMel133 cells, however, the cell viability IC50s for RAFi and MEKi are at least one order of magnitude higher than the doses required to reduce the phosphorylation of immediate downstream targets (phosphoMEK and phosphoERK) by 50% (Figure 6B, Table 2), and hence, the cells are considered to be relatively resistant to both drugs. When combined with JQ1, the concentrations of RAFi/MEKi that are sufficient to reduce cell viability by 50% are in the range of MEK and ERK phosphorylation IC50s (Table 2). This is a direct result of the strong and synergistic increase in G1arrest phenotype in response to the drug combinations (Figure 6D). Thus, the JQ1MEKi/RAFi combinations shift the required doses to induce phenotypic responses close to the doses required for inhibition of the ERK pathway activity. Interestingly, in an independent study in hematopoietic cancer, the combination of JQ1 affecting the cMyc levels and a FLT3 tyrosine kinase inhibitor (TKI) was synergistic in overcoming resistance to the FLT3 TKI in acute myeloid leukemia (AML) cells (Fiskus et al., 2014). As far as we know, our study is the first demonstration of a combined inhibition of a BET bromodomain (BETBRD) protein and a protein kinase molecule to overcome drug resistance in solid tumors (Filippakopoulos and Knapp, 2014).
To interpret these positive results, we note that cMyc is a multifunctional transcription factor with a central role in cell cycle progression and receives inputs from diverse signaling pathways, including ERK, PI3K/AKT, Wnt, and TGFβ signaling (Dang, 2012). We hypothesize that, beyond the efficacy of cotargeting BRAF kinase and cMyc, there may be a particular advantage in blocking proliferative signals from various oncogenic pathways that converge on the downstream effector cMyc. In other words, cMyc may be a particular ‘point of vulnerability’ in RAFi drugresistant cells (Figure 7) and this may be reflected in the high rank of cMyc as a cotarget to induce G1arrest phenotype in the datadriven perturbation biology network models.
Although the experiments here are limited to the SkMel133 melanoma cell line, we are tempted to put forward a more general hypothesis. We suggest that synergistic cotargeting of a key genomic aberration (e.g., BRAFV600E) and a downstream point of convergence of signaling (e.g., cMyc) may be a more effective way of overcoming the resistance than cotargeting another upstream pathway (Figure 7). While inhibition of the oncogenic activation by the first drug can be compensated by alteration in a number of bypass pathways, it is plausible that most if not all signaling events have to converge on similar, proproliferative consequences—in this example via the downstream transcriptional activator cMyc. Inhibition of a downstream signaling event may therefore be a generally more effective cotargeting strategy for an entire set of alternative upstream bypass mechanisms of resistance. The key proviso, of course, is absence of toxicity, which may be more likely for cotargeting the upstream oncogenic activator (not active in normal cells, e.g., BRAV600E) in combination with a downstream effector (functional in both normal and tumor cells, e.g., cMyc [Soucek et al., 2008]) at optimized drug concentrations. In other words, the upstream–downstream combination of targets may have the dual advantage of generality in overcoming resistance and avoidance of toxicity. To generally test this hypothesis, one would have to go beyond the scope of this work and perform cotargeting experiments in other melanoma cell lines, for which the mechanism of resistance is known, or more generally, in other cell lines or mouse models of resistance to targeted therapeutics.
In the current study, we have only tested the JQ1 and RAFinhibitor combinations in vitro. Future in vivo studies will be important for identifying the optimal BETBRD inhibitor for combination therapy in melanoma and other tumors. This may be clinically relevant, as BETBRD inhibitors hold substantial clinical potential for treatment of both solid and hematological tumors. Indeed, the observation that BETBRD inhibitors (e.g., JQ1) effectively block cellular growth through transcriptional repression of cMyc in mouse models has created the rationale for clinical trials in cMycdependent cancer types (Mertz et al., 2011). Unfortunately, the widely used investigational BETBRD inhibitor JQ1 has been reported to have poor in vivo efficacy with low metabolic stability and biosolubility, but fortunately, other BETBRD inhibitors with more promising pharmacokinetic properties have been and are being developed. Specifically, JQ1 is reported to have a maximum blood concentration of 0.84 μM in mouse models, whereas the compounds IBET762 and IBET151 reach maximum blood concentrations of 26 μM and 82 μM, respectively (Dawson et al., 2011). In general, triazolobenzodiazepine (e.g., IBET762 and IBET151 compounds) and triazolothienodiazepine (e.g., OTX015, TEN010, CPI0610 compounds) BETBRD inhibitors are reported to be stable in liver microsomes and have up to 30% bioavailability. These compounds may therefore be more suitable for in vivo applications (Filippakopoulos and Knapp, 2014). Several clinical trials of BETBRD inhibitors with desirable pharmacokinetic properties have been initiated in multiple cancer types, including trials with OTX015 for the treatment of acute leukemia (ClinicalTrials.gov identifier: NCT01713582), CPI0610 for progressive lymphoma (ClinicalTrials.gov identifier: NCT01949883), IBET762 and TEN010 for NUT midline carcinoma (ClinicalTrials.gov identifier: NCT01587703, NCT01987362), and a dose escalation and cohort expansion study for TEN010 in patients with AML and myelodysplastic syndrome (ClinicalTrials.gov identifier: NCT02308761). The trials with BETBRD inhibitors in different cancer subtypes suggest use of such compounds either as single agents or in combinations may be a plausible strategy for treatment of other cancer types, including RAFinhibitor resistant melanoma.
The modelbased predictions provide comprehensive and testable hypotheses on complex signaling mechanisms relevant to the development of novel therapies. More comprehensive experimental data and improvements in signaling databases are likely to lead to network models with even higher predictive power. Interpreting the cell linespecific predictions on the background of genomic profiling of tumor samples from patients may guide therapeutic choices for tumors with similar genomic backgrounds. To assist in these efforts, we have developed genomics methods to classify tumors based on patterns of functional oncogenic alterations and to identify cell lines that most closely resemble tumors (Ciriello et al., 2013; Domcke, Sinha, et al., 2013). By integrating network models, genomics, and pathway analysis, one can expect to generate wholecell models of signaling and drug response in mammalian cells with potential applications to personalized medicine and genomically informed clinical trials.
Website
Website for the perturbation biology method: http://www.sanderlab.org/pertbio/. The website stores the source code for BPdecimation algorithm, input scripts, and data files for running the BPguided decimation code, perturbation response data, downloadable files for the executable model solutions, and the complete simulation results. A detailed description of the perturbation biology method is also provided in the webpage.
Materials and methods
Computational methods
Automated extraction of prior information from signaling databases
Request a detailed protocolPathway Extraction and Reduction Algorithm (PERA) was developed to automatically extract prior information from multiple signaling databases and generate a prior information network. The input to PERA is a list of (phospho) proteins identified by their HGNC symbols (e.g., AKT1), phosphorylation sites (e.g., pS473), and their molecular status (i.e., activating or inhibitory phosphorylation, total concentration). The output of PERA is a set of directed interactions between signaling molecules represented in a Simple Interaction Format (SIF). The PERA software is available at http://bit.ly/bp_prior as a free software under LGPL 3.0.
Requirements for prior information
There is a growing number of signaling databases (Bader et al., 2006) that capture pathway information in high detail. Compared to the interaction networks with relatively low level of detail, pathway databases contain information about the phosphorylation states of proteins and interactions specific to these states. Although this corpus is highly valuable for phosphoproteomic analysis, it remains mostly untapped. To utilize this knowledge, one needs algorithms to map the experimentally measured proteins and phosphoproteins to their pathway database equivalents, find connections between them, and reduce the outcome to a format that can be used for inference and analysis.
The algorithm
We developed the algorithm and the software tool, PERA to automatically extract prior information from multiple signaling databases in the BioPAX (Demir et al., 2010) format and generate a prior information network. PERA takes a list of (phospho)proteins identified by their HGNC symbols (e.g., AKT1), phosphorylation sites (e.g., pS473), and their molecular status (i.e., activating or inhibitory phosphorylation, total concentration) as input and then finds directed signaling paths between these entities. These paths are then reduced to directed interactions between signaling molecules represented in a SIF.
The PERA framework
The prior information is extracted from components of the Pathway Commons 2 database in three steps: first, using the pathsbetween graph query algorithm (Dogrusoz et al., 2009), PERA generates a subgraph (i.e., graphofinterest) of Pathway Commons, which contains all the input proteins and all known connections within their first neighborhoods. Second, using the phosphorylation and activity state information, input entities are mapped to the corresponding protein states in the graphofinterest. During this mapping step, protein states that do not match with either the corresponding annotation for phosphorylation or activity state are filtered out. Phosphorylation site mismatches up to six residues are tolerated during the filtering step to account for phosphorylation site ambiguities due to either database curation errors or crossorganism annotations. Third, paths that result in the addition or removal of a phophogroup at a phosphorylation site are extracted and mapped to phosphoprotein nodes. For total protein nodes, all nonphosphoprotein specific, directed signaling paths are included. For this application, the maximum allowable graph query distance limit was set to 1 and only the Reactome (Matthews et al., 2009) and NCINature PID pathway (Schaefer et al., 2009) data resources were used. Although we limited ourselves to short path distances and two pathway databases that we were most familiar with, PERA can be applied to extract information from any pathway database that exports to BioPAX and can be configured for searching paths of arbitrary length.
Key considerations for extracting pathway information
Two major semantic issues need to be taken into account when extracting statespecific pathway information from databases. First semantic issue stems from the ambiguities in mapping proteins with multiple phosphorylation sites. When a protein with multiple phosphorylation sites is experimentally profiled with an antibody, which recognizes only a single phosphorylation site, the antibody will actually bind to a heterogeneous mixture of phosphostates provided that the epitope is phosphorylated (e.g., antiAKTpS473 Ab may bind to both AKTpS473 and AKTpS473_pT308 but not to AKTpT308). For proteins with multiple observed phosphorylation sites, this might lead to semantic conflicts since a double phosphorylated node should be mapped to both observations (i.e., single and doublephosphorylated states). We included an optional ‘strict’ mapping scheme to map only the phosphoproteins that exactly match the observed node—in our case always single phosphoproteins (e.g., the epitope of antiAKTpS473 Ab is mapped to AKTpS473 but not to AKTpS473_pT308). Since our extraction algorithm is much more tolerant of missing interactions compared to false interactions, we opted to use this flag for this application.
The second semantic issue stems from the fact that pathway databases are often curated from multiple independent studies, spanning multiple cellular states, cell and tissue types, and even multiple model organisms. As a result, they are a superimposition of possible interactions over a wide range of spatiotemporal and genetic contexts. On the other hand, databases cover only a subset of all possible contexts. In our case, we expect only a subset of the interactions in the pathway databases to be active in our cell lines and cover only a subset of the interactions that we observe. This observation necessitates incorporating prior information as ‘soft’ restraints for network inference. For this purpose, we devised a modified cost function, which includes a term for prior information (See below. Also see [Molinelli, Korkut, Wang et al., 2013] and [Miller et al., 2013]).
Mathematical description of network models
Request a detailed protocolThe network models represent the time behavior of the cellular system in a set of perturbation conditions as a series of coupled nonlinear ODEs (Nelander et al., 2008).
Equation 1. Network model ODEs
In the network models, each node represents the quantitative change of a biological variable, x^{μ}_{i} ((phospho)protein level and phenotypic change) in the perturbed condition, μ, relative to the unperturbed condition. W_{ij} quantifies the edge strength, which is the impact of upstream node j on the time derivative of downstream node i. We assign semidiscrete values to each W_{ij}, W = {w_{ij}, ∀w_{ij} ∈ {−1, −0.8,…,0.8, 1}}. α_{i} constant is the tendency of the system to return to the initial state, and ε_{i} constant defines the dynamic range of each variable i. The transfer function, Φ(x) ensures that each variable has a sigmoidal temporal behavior.
Modified cost function for network inference with prior information
Request a detailed protocolWe have quantified the cost of a model solution by an objective cost function C(W). The network configurations with low cost represent the experimental data more accurately. Here, we have incorporated an additional prior information term to the cost function to construct models with improved predictive power. The newly introduced term in the cost function accounts for the prize introduced when the inferred w_{ij} is consistent with the prior information.
In order to generate predictive network models of signaling, we have quantified the cost function C(W) to penalize (a) discrepancies between predicted (x_{i}^{μ}) and experimental (x_{i}^{μ*}) values of the system variables at a time point {t_{l}} in condition μ, and (b) the number of nonzero interactions with an L0 norm (Nelander et al., 2008; Molinelli, Korkut, Wang et al., 2013). Here, we have also incorporated prior information to construct network models with higher predictive power even if the models are constructed with limited experimental data. The newly introduced third term in the cost function accounts for the prize introduced when the inferred w_{ij} conforms to the prior information (w_{ij}^{prior}). The modified cost function with prior information term is formulated as in Equation 2.
Equation 2. Modified cost function
In Equation 2, the first term accounts for the discrepancies between predicted (x_{i}^{μ}) and experimental (x_{i}^{μ*}) values of the system variables at a time point {t_{l}} in condition μ, for a particular network configuration (W = {w_{ij}}). The second term is the complexity factor with an L0 norm, which reduces the number of nonzero interactions in a network configuration and ensures that resulting network models are sparse. The final term, η(W) is the cumulative prior prize function for W = {w_{ij}}. η(w_{ij} = ω) has a negative real value and reduces the overall model cost if the inferred ω conforms to the prior information. Here, we will formulate the newly introduced prior information term in the modified cost function.
Generalized prior information prize function
In order to derive the prior information prize (PIP) function, we first assume the value of prior w_{ij} is normally distributed (truncated by boundaries of w_{ij}) around an expected prior value, E[w_{ij}^{prior}] with a standard deviation, σ_{ij}^{prior}. The prior probability model for each interaction in the prior network is formulized using the normal distribution probability density function around E[w_{ij}^{prior}] as in Equation 3.
Equation 3. Probability model of prior observations
The Equation 3 implies that the probability of inferring a particular w_{ij} value and consequently, the prize introduced to the cost function increase as w_{ij} approaches E[w_{ij}^{prior}]. We calculate the prize for the fit between the inferred w_{ij} and the E[w_{ij}^{prior}] with an inverse Boltzmann conversion of the probability model (Jaynes, 1957).
Equation 4. Error model of individual prior observations
The first term introduces the maximum possible prize for each interaction that is represented in the prior information network. Second term penalizes the discrepancies between the w_{ij} and the E[w_{ij}^{prior}]. The penalty for the discrepancy is zero when w_{ij} = E[w_{ij}^{prior}] and the prior prize, η(w_{ij}) assumes the highest possible value. ${\mathrm{\kappa}}_{\text{ij}}\left(\text{ln}\left(1/\left({\mathrm{\sigma}}_{\text{ij}}^{\text{prior}}\sqrt{2\pi}\right)\right)\right)$ is a constant analogous to the inverse temperature in statistical physics. The constant κ_{ij} and the resulting η(w_{ij}) are zero if the interaction between nodes i and j is not included in the prior network model. By summing over all i and j, the cumulative prior prize for all W = {w_{ij}} becomes Equation 5.
Equation 5. Cumulative error model of prior observations
Even though the current PIP function is derived with the normal distribution assumption, a variety of alternative distributions such as a bimodal distribution can be fitted in Equation 3. Such efforts will lead to a variety of customdesigned PIP functions when necessary.
Estimation of E[w_{ij}^{prior}] and σ_{ij}^{prior}
An important requirement for the implementation of Equation 5 is the accurate estimation of E[w_{ij}^{prior}] and σ_{ij}^{prior}. Here, we propose two estimation strategies that are alternative to systematic extraction of priors from databases.
Priors from experiments
Our first strategy relies on use of biochemical measurements with multiple biological replicates to quantify the influence of node j on node i (w_{ij}). The E(w_{ij}^{prior}) and σ_{ij}^{prior} can be estimated from the mean value (<w_{ij}^{prior}>) and standard deviation (σ_{ij}) of the readouts from biological replicates in such measurements. Depending on the nature of the interaction under study, a variety of different biochemical methods can be utilized. These methods include but are not limited to enzyme activity assays, quantification of protein–protein interactions (Selbach and Mann, 2006), or additional RPPA readouts under specific perturbation conditions.
Priors from network models
As we build network models of signaling in diverse biological contexts using perturbation data and inference algorithms, a vast amount of information on signaling properties will become available. Consequently, we will be able to estimate the E(w_{ij}^{prior}) and σ_{ij}^{prior} for newly inferred models from available network models. Previously inferred {w_{ij}} values and standard deviation in BPgenerated probability distributions will serve as a basis to estimate the E(w_{ij}^{prior}) and σ_{ij}^{prior}, respectively.
The simplified PIP function
The current form of the prior information is limited to a set of binary interactions due to the qualitative nature of the databases, from which we extract the information. Therefore, we used a simplified PIP function to introduce the prior information restraints to the inference scheme. The prior information stored in the databases may correspond to activating (w_{ij} > 0), inhibitory (w_{ij} < 0), or generic (w_{ij} ≠ 0) interactions. In such situations, w_{ij}^{prior} has a uniform distribution within the defined boundary. Thus, the E[w_{ij}^{prior}] can be set as equal to any value within the interval, [w_{ij}^{minprior}, w_{ij}^{maxprior}]. For example, a prior for a positive interaction between nodes i and j is represented by an interval of (0, w_{ij}^{max}), where w_{ij}^{max} is the maximum allowed edge strength in the inference scheme. Consequently, the squared distance error in Equation 4 vanishes and the prior prize assumes a fixed value over the defined prior interval.
Equation 6. Simplified error model of prior information
The binary nature of the interactions also implies a generic κ value $\left(\mathrm{\kappa =}{\mathrm{\kappa}}_{\text{ij}}\left(\text{ln}\left(1/\left({\mathrm{\sigma}}_{\text{ij}}^{\text{prior}}\sqrt{2\pi}\right)\right)\right){,\text{\hspace{0.17em}\kappa}}_{\text{ij}}>0,\forall \text{\hspace{0.17em}i},\text{j}\in \text{N}\right)$ for all interactions represented in prior information network and κ = 0 for all other situations. The resulting prior information term in the cost function is a step function as shown in the main text methods section.
Equation 7. Simplified error model of individual prior observations
For each interaction in the prior information network, η(w_{ij}) has a negative value. Thus, the prior information introduces a prize in the overall cost function. The formulation and implementation into the BP equations create a softrestraint for prior information extracted from signaling databases. In our probabilistic inference framework, a prior information is accepted as an edge (wij ≠ 0) only if the prize introduced can override any discrepancy between predicted (x_{i}^{μ}) and experimental (x_{i}^{μ*}) values of the system variables. Thus, the contextspecific character of inferred models is preserved when the cost function in Equation 2 is used.
Network model construction and response prediction
Request a detailed protocolNetwork models are constructed with a twostep strategy. The method is based on first calculating probability distributions for each possible interaction at steady state with the BP algorithm and then computing distinct solutions by sampling the probability distributions. We described the theoretical formulation, the underlying assumptions and simplification steps of the BP algorithm for inferring network models of signaling elsewhere (Molinelli, Korkut, Wang et al., 2013). The network models include 82 proteomic, 5 phenotypic, and 12 activity nodes. Activity nodes couple the effect of drug perturbations to the overall network models (See Molinelli, Korkut, Wang et al. for quantification of activity nodes).
Belief propagation
Belief propagation (BP) algorithm iteratively approximates the probability distributions of individual parameters. The iterative algorithm is initiated with a set of random probability distributions. In each iteration step, individual model parameters are updated (local updates) based on the approximate knowledge of other parameters, experimental constraints, and prior information (global information). In the next iteration, the updated local information becomes part of the global information and another local update is executed on a different model parameter. The successive iterations continue over different individual parameters until the updated probability distributions converge to stable distributions. The iterations between the local updates and the global information create an optimization scheme that W = {w_{ij}} is inferred given a probability model. Explicitly, the following cavity update equations are iteratively calculated until convergence.
Equation 8. BP update equations
In Equation 8a, P^{μ}(w_{ij}) approximates the mean field of the parameters with a sparsity constraint (λδ(w_{ij})) and a bias from prior information restraints (η(w_{ij})). In Equation 8b, ρ^{μ}(w_{ij} = ω) is a mean field derived change to the probability distribution of the locally optimized parameter, towards minimizing the model error (C^{SSE}(W)). (See Molinelli, Korkut, Wang et al. (2013) for derivation and implementation of BP equations).
Gaussian integration of cavity parameter update
Information from a sufficiently large set of noncavity constraints and parameters is integrated to update the cavity parameters. As a corollary to the central limit theorem, we can assume that the meanfield effect of noncavity parameters on the cavity update follows a Gaussian distribution. The mean and variance of the Gaussian distribution are the mean values and variances of the global distributions given by individual distributions of noncavity P^{μ}(w_{ij}) ∀j ≠ k. Thus, we replace the sum over multivariate configurations of all noncavity parameters with a single Gaussian integration over the meanfield of noncavity configurations.
Equation 9. Gaussian approximation to the cavity parameter update
See Equations 4–14 in Molinelli, Korkut, Wang et al. (2013) for derivation of Equation 9 and detailed descriptions of the terms in the equation. Note that, BP parameters (temperature factor, β = 2; complexity factor; λ = 5; prior information weight, κ = 5) are adapted for all computations after a systematical analysis of all three parameters (Figure 4—figure supplement 5). When prior information is used to construct the models, the complexity term, λ is set to 2.5 for phenotypic nodes in order to approximately match the complexity around the phenotypic nodes and rest of the networks. When the parameter space is sufficiently large, the Gaussian approximation significantly reduces the computational complexity of the problem without loss of overall accuracy. By means of the iterative approach with the meanfield Gaussian approximation, the timecomplexity of the problem is strongly reduced and the obstacles imposed by combinatorial complexity are circumvented.
BPguided decimation
Distinct networks models are instantiated from BPgenerated probability distributions with the BPguided decimation algorithm (Figure 1—figure supplement 1) (Montanari et al., 2007). This procedure generates distinct and executable network models. In this study, 4000 distinct network models are generated in each computation.
Simulations with in silico perturbations
Network models are executed with specific in silico perturbations until all system variables {x_{i}} reach steady state. The perturbations acting on node, i, are exerted as realvalued, u_{i}^{μ}, vectors in model Equation 1. The DLSODE integration method (ODEPACK) (Hindmarsh, 1993) is used in simulations (default settings with, MF = 10, ATOL = 1e10, RTOL = 1e20).
Experimental methods
Cell cultures and perturbation experiments
Request a detailed protocolAll the perturbation experiments are performed using the RAFinhibitor resistant SkMel133 cell line. The identity of the SkMel133 was checked for mislabeling, contamination, and misidentification using a multiplexed PCR/mass spectrometry–based genetic fingerprinting assay run on the HapMap Sequenom platform. In perturbation experiments, the drugs are applied to cell cultures after SkMel133 cells are grown to about 40% confluence in complete RPMI1640 medium (10% heatinactivated fetal bovine serum, 100 units/ml each of penicillin and streptomycin, and incubated at 37°C in 5% CO_{2}) in 6well plates. 24 hr after drug administration, the perturbed cells are harvested. In control experiments (i.e., no drug condition), cells are treated with the DMSO drug vehicle for 24 hr. SkMel133 cells are perturbed with 12 targeted drugs applied as single agents or in paired combinations (see Supplementary file 1A,C for drug list, presumed targets, dosing and sources). In total, cells are treated with 89 unique perturbations. In paired combinations, each drug concentration is selected to inhibit the readout for the presumed target or the downstream effectors by 40% (IC40) as determined by Western blot experiments (Molinelli, Korkut, Wang et al., 2013) (Supplementary file 1A). In single agent perturbations, each drug is applied at two different concentrations, IC40 and 2 × IC40. In validation experiments, (+)JQ1 (Cayman Chemicals, Ann Arbor, MI) and the U.S. Food and Drug Administration (FDA)approved RAFi PLX4032 (Selleckchem, Houston, TX) are used.
Reverse phase protein arrays
Request a detailed protocolProteomic response profiles to perturbations are measured using RPPA (MD Anderson Cancer Center RPPA Core Facility) (Tibes et al., 2006). The cells are lysed 24 hr after drug treatment. For sample preparation, cell pellets are lysed in RPPA buffer (1% Triton X100, 50 mM HEPES pH7.4, 150 mM NaCl, 1.5 mM MgCl_{2}, 1 mM EGTA, 100 mM NaF, 10 mM Na_{4}P_{2}O_{7}, 1 mM Na_{3}VO_{4}, 10% glycerol, freshly added protease and phosphatase inhibitors). After cells are lysed, total protein concentrations in cell lysates are measured with the Bradford assay and the final protein concentrations are adjusted to 1–1.5 mg/ml. Samples are boiled in SDS sample buffer without bromophenol blue for 5 min at 95°C. Cell lysates are spotted on nitrocellulosecoated slides as described in Tibes et al. (2006). Antibody staining intensities are quantified using the MicroVigene automated RPPA module (VigeneTech, Inc., Carlisle, MA) and the standard RPPA protein concentration normalization procedure (Neeley et al., 2009) is followed. Three biological replicates are spotted for each sample (i.e., drug condition) on RPPA slides. Each slide is interrogated with a particular antibody, so for each experimental condition 138 proteomic entities (levels of total protein or phoshoprotein) are profiled on 138 slides (Supplementary file 1C). The proteomic readouts are log normalized with respect to the corresponding untreated condition readouts. We have eliminated those readouts with intra or interslide coefficient of variation >0.15 (i.e., low reproducibility) and low degree of staining by antibodies. 100 proteomic entities are chosen for further analysis.
Detection of proteomic entities that respond to perturbations
Request a detailed protocolTo generate a constraint set for network inference, we excluded the proteomic entities that do not respond to at least a single perturbation condition from the network models. An iterative signaltonoise detection algorithm is used for this purpose. At each iteration, the standard deviation of the overall data is computed and data points that deviate for >2.5σ_{i} (σ_{i}: standard deviation of all data points at iteration i) from the data average are selected as signal (i.e., responding to perturbation). In the next iteration, the σ_{i + 1} is computed for all data points except the data points selected as signal in the previous iteration. The σ calculation is iterated until convergence. The proteomic entities that respond to at least a single perturbation condition are selected. This has lead to highly reliable response measurements on a final set of 82 proteomic entities, which have been used in the network models.
Quantitative phenotypic assays
Request a detailed protocolAll phenotypic measurements are made in perturbation conditions identical to those of proteomic measurements. Cell viability and cell cycle progression are measured using the resazurin assay (72 hr after drug treatment) and flow cytometry analysis (24 hr after drug treatment), respectively. The percentage of cells in the G1, G2/M, and S phases and subG1 fraction are recorded based on the respective distribution of DNA content in each phase. For the resazurin assay, each well is treated with the resazurin (Sigma–Aldrich, St. Louis, MO, Catalog #R7017) at a final concentration of 44 μM for 1 hr. The fluorescent signals are measured at 530nm excitation wavelength and 590nm emission wavelength. Standard curves of cell numbers are used to calculate the cell numbers in each sample. The cell cycle progression phenotypes are measured with flow cytometry analysis. Cells are seeded in 6well plates. 72 hr after drug application, cells are harvested by trypsinization, including both suspended and adherent fractions, and washed in cold PBS. Cell nuclei are prepared as described in Nusse and Kramer (1984), and cell cycle distribution was determined by flow cytometric analysis of DNA content using red fluorescence of 488nm excited ethidium bromidestained nuclei. The percentage of cells in the G1, G2/M, and S phases and subG1 fraction are recorded.
Advances in the perturbation biology method
The improvements in perturbation biology have enabled us to establish our method as a widely applicable tool in genomically informed preclinical studies. Here, we obtained signaling models with increased scope (increased coverage of diverse pathways using a data set 20fold bigger than previously used data sets) and nominated a particular drug combination to overcome the drug resistance in melanoma. The increased scope of the models is particularly critical for biologists to adopt our method as such models provide a nonreductionist insight into relations between the molecular changes and complex biological phenomena such as drug response, cell fate in diverse conditions, or any other cellular phenotype that can be measured quantitatively. Here, we improved the perturbation biology method and solved three key challenges in network inference to reach scopes and predictive power necessary for addressing complex problems in cancer biology.
A biologically grounded modeling framework that can be incrementally improved
Request a detailed protocolOur current modeling scheme can incorporate both quantitative and qualitative prior information through the introduction of a generalized probabilistic error model to the inference cost function (Equations 2–7). In the generalized form of the error model, the prize introduced by the priors is sampled from a probability model with a Gaussian distribution (Equation 3) for the strength (w_{ij}^{prior}) of each prior interaction. A Boltzmann conversion of the prior probability model leads to an error model that is fully compatible with the BPgenerated probability distributions, P(w_{ij}). Such a generalized error model will be most useful when quantitative prior information is available for the signaling interactions.
The general error model provides an advantage in building celltype specific models in diverse but related biological samples. When we build network models for a series of biological systems (e.g., cell lines or patient derived samples) with similar but distinct genomic properties, the BPgenerated probability distributions in one cellular system can be extrapolated to another system as a set of probabilistic biases (i.e., priors). The probabilistic extrapolation between different biological samples will significantly reduce the amount of experimental constraints (i.e., perturbation conditions) required to build predictive models while preserving the celltype specificity for each biological sample studied. This is extremely promising in the context of cancer biology where a model portfolio built from a cohort of patientderived samples (cell lines or mouse avatars) coupled with genomic profiling can inform treatment decisions for new patients. This bridge between genomic profiles, models, and response to treatment is an open challenge that we think is going to be critical for future precision medicine applications.
The versatile prior error model introduced here enables us to use known biology in various forms that can be available from previous studies or databases. On one extreme, the error model is simplified and reduced to a step function (Equation 7). The simplified form of the error model is compatible with the binary priors as are the priors used in this study, which were extracted from the curated signaling databases. Alternatively, one can generate hybrid forms of the cumulative prior error model. Such hybrid models may include (i) the truncated Gaussian model for a subset of priors with quantitative parameters, (ii) the simplified step function form for binary priors, and (iii) alternative forms (e.g., based on a bimodal distribution) for other priors. Although we used the simplified form of the error model and constructed models in a single cell line, we believe the flexibility of the prior error model will be transformative for generating highly predictive models of different biological systems. Such models will subsequently capture the intersample variations in response to combinatorial perturbations.
Systematic integration of prior information
Request a detailed protocolDevelopment of the PERA algorithm leads to systematic extraction and use of prior interactions, which has a series of advantages over ad hoc introduction of priors in the network modeling scheme. First, currently available network inference algorithms and prior network extraction utilities work at either gene or protein level. Therefore, they are of limited use for inference algorithms that operate on phosphoprotein levels. PERA allows users to specify both phospho and totalprotein levels of the genes and provide users with prior knowledge with higher granularity. Second, pathway resources are regularly being updated with newly curated interactions; therefore, our prior knowledge on entities of interest is constantly evolving. PERA can work against new versions of pathway data from Pathway Commons and this enables users to obtain uptodate prior information for their proteins by simply rerunning the application. Third, PERA extracts prior information from pathway data integrated from multiple resources in an automated manner; hence, the interactions extracted from the databases are unbiased compared to the prior knowledge obtained from literature by a single researcher. The algorithm also allows users to specify options (such as the level of phosphorylation site mismatch tolerance and the maximum length between two entities to be considered) for adjusting the complexity and the specificity of the final prior network. Here, an open challenge for future is assigning quantitative confidence scores and statistical information for each prior interaction extracted from databases. Currently, the databasederived priors are in binary form. However, with incorporation of literature cocitation data and experimental coexpression data in different tissues will enable us to estimate the confidence scores and statistical parameters for potential prior interactions.
Performance improvements that can scale the method to our experimental capabilities
Request a detailed protocolWe optimized the performance of our algorithm by introducing heuristics at key bottlenecks. This was necessitated by 20fold increase in data points compared to our previous work. One of the key bottlenecks in the modeling scheme is the BPguided decimation, for which the BPderived probability distributions are updated after each step (i.e., after sampling all the BPderived probability distributions and fixing the interaction value for one edge with a probabilistic draw) (see supplementary Figure 1—figure supplement 1). In other words, the posterior probability of a W_{ij} is converted to 1. Not surprisingly, the BPprobability for significant fraction of the selected edges is within the range of 0.95–0.99. As fixing an individual edge from P^{unfixed}(W_{ij}) = 0.95 (P^{unfixed}: probability before fixing the edge value) to P^{fixed}(W_{ij}) = 1 does not significantly alter the probabilistic landscape of the global solution space that involves thousands of other probabilistic events, we updated the global interaction landscape after each decimation step if the P^{unfixed}(W_{ij}) < 0.95 for the selected edge. However, accumulation of such residual shifts in individual probabilities may lead to a significant change in the global landscape after consecutive iterations of the decimation algorithm. To prevent such inaccuracy in modeling, we updated the probability landscape after five consecutive decimation steps even if P^{unfixed}(W_{ij}) > 0.95 at each step. In practice, the conditional updating scheme described here provided a >threefold increase in the performance of the BPdecimation algorithm without substantially modifying the sampled probability space.
References

1
WNT5A enhances resistance of melanoma cells to targeted BRAF inhibitorsThe Journal of Clinical Investigation 124:2877–2890.https://doi.org/10.1172/JCI70156

2
Pathguide: a pathway resource listNucleic Acids Research 34:D504–D506.https://doi.org/10.1093/nar/gkj126
 3
 4
 5

6
Pathway commons, a web resource for biological pathway dataNucleic Acids Research 39:D685–D690.https://doi.org/10.1093/nar/gkq1039

7
Improved survival with vemurafenib in melanoma with BRAF V600E mutationThe New England Journal of Medicine 364:2507–2516.https://doi.org/10.1056/NEJMoa1103782
 8
 9

10
Emerging landscape of oncogenic signatures across human cancersNature Genetics 45:1127–1133.https://doi.org/10.1038/ng.2762

11
ARN509: a novel antiandrogen for prostate cancer treatmentCancer Research 72:1494–1503.https://doi.org/10.1158/00085472.CAN113948
 12
 13
 14
 15

16
The BioPAX community standard for pathway data sharingNature Biotechnology 28:935–942.https://doi.org/10.1038/nbt.1666
 17

18
Evaluating cell lines as tumour models by comparison of genomic profilesNature Communications 4:2126.https://doi.org/10.1038/ncomms3126

19
Targeting bromodomains: epigenetic readers of lysine acetylationNature Reviews. Drug Discovery 13:337–356.https://doi.org/10.1038/nrd4286
 20

21
Systems biology and combination therapy in the quest for clinical efficacyNature Chemical Biology 2:458–466.https://doi.org/10.1038/nchembio817

22
Combined BRAF and MEK inhibition in melanoma with BRAF V600 mutationsNew England Journal of Medicine 367:1694–1703.https://doi.org/10.1056/NEJMoa1210093

23
Improved survival with MEK inhibition in BRAFmutated melanomaThe New England Journal of Medicine 367:107–114.https://doi.org/10.1056/NEJMoa1203421
 24

25
Genetic alterations in signaling pathways in melanomaClinical Cancer Research 12:2301s–2307s.https://doi.org/10.1158/10780432.CCR052518

26
ODEPACK, a systematized collection of ODE solversIn: RS Stepleman, M Carver, R Peskin, WF Ames, R Vichnevetsky, editors. Scientific Computing. Amsterdam: NorthHolland Publishing Company. pp. 55–64.
 27

28
Information theory and statistical mechanicsPhysical Review 106:620–630.https://doi.org/10.1103/PhysRev.106.620
 29

30
Network quantification of EGFR signaling unveils potential for targeted combination therapyMolecular Systems Biology 9:673.https://doi.org/10.1038/msb.2013.29

31
The systems biology graphical notationNature Biotechnology 27:735–741.https://doi.org/10.1038/nbt.1558
 32

33
Sensitivity of human lung adenocarcinoma cell lines to targeted inhibition of BET epigenetic signaling proteinsProceedings of the National Academy of Sciences of USA 109:19408–19413.https://doi.org/10.1073/pnas.1216363109

34
Combined BRAF and MEK inhibition versus BRAF inhibition alone in melanomaThe New England Journal of Medicine 371:1877–1888.https://doi.org/10.1056/NEJMoa1406037
 35

36
Reactome knowledgebase of human biological pathways and processesNucleic Acids Research 37:D619–D622.https://doi.org/10.1093/nar/gkn863

37
Targeting MYC dependence in cancer by inhibiting BET bromodomainsProceedings of the National Academy of Sciences of USA 108:16669–16674.https://doi.org/10.1073/pnas.1108190108
 38

39
Perturbation biology: inferring signaling networks in cellular systemsPLOS Computational Biology 9:e1003290.https://doi.org/10.1371/journal.pcbi.1003290

40
Solving constraint satisfaction problems through belief propagationguided decimationarXiv, 0709.1667. http://arxiv.org/abs/0709.1667.
 41

42
Variable slope normalization of reverse phase protein arraysBioinformatics 25:1384–1389.https://doi.org/10.1093/bioinformatics/btp174

43
Models from experiments: combinatorial drug perturbations of cancer cellsMolecular Systems Biology 4:216.https://doi.org/10.1038/msb.2008.53
 44
 45

46
Causal network inference using biochemical kineticsBioinformatics 30:I468–I474.https://doi.org/10.1093/bioinformatics/btu452
 47
 48
 49

50
(V600E)BRAF is associated with disabled feedback inhibition of RAFMEK signaling and elevated transcriptional output of the pathwayProceedings of the National Academy of Sciences of USA 106:4519–4524.https://doi.org/10.1073/pnas.0900780106

51
Targeting MYCN in neuroblastoma by BET bromodomain inhibitionCancer Discovery 3:308–323.https://doi.org/10.1158/21598290.CD120418

52
Molecular markers of tumor progression in melanomaCurrent genomics 10:231–239.https://doi.org/10.2174/138920209788488526
 53
 54

55
PID: the pathway interaction databaseNucleic Acids Research 37:D674–D679.https://doi.org/10.1093/nar/gkn653

56
A shrinkage approach to largescale covariance matrix estimation and implications for functional genomicsStatistical Applications in Genetics and Molecular Biology 4:Article32.https://doi.org/10.2202/15446115.1175
 57
 58
 59

60
Survival in BRAF V600mutant advanced melanoma treated with vemurafenibThe New England Journal of Medicine 366:707–714.https://doi.org/10.1056/NEJMoa1112302
 61
 62
 63
 64
 65
 66
 67

68
Intratumor heterogeneity alters most effective drugs in designed combinationsProceedings of the National Academy of Sciences of USA 111:10773–10778.https://doi.org/10.1073/pnas.1323934111
Decision letter

Edoardo M AiroldiReviewing Editor; Harvard University, United States
eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.
Thank you for sending your work entitled “Perturbation biology models predict cMyc as an effective cotarget in RAF inhibitor resistant melanoma cells” for consideration at eLife. Your article has been favorably evaluated by Aviv Regev (Senior editor), a Reviewing Editor, and two additional reviewers.
The Reviewing editor and the reviewers discussed their comments before we reached this decision, and the Reviewing editor has assembled the following comments to help you prepare a revised submission.
The editorial team has identified two areas of major concern: (1) novelty of the approach and (2) biological significance.
Regarding novelty of the approach, a portion of the most interesting methodological innovations has been already published by a subset of the authors. The preprocessing PERA step is new, and so are a number of other technical advances. In this paper, the authors extended the scope of the analysis substantially (phosphoprotein antibodies from 16 to 138, perturbation conditions from 44 to 89, and phenotypes from 1 to 5), and we appreciate that this extension is critical for biologists to adopt such approaches. However, the authors fail to clearly present the innovations over previous work. The authors should make an effort, perhaps in a dedicated subsection under Methods, to offer a nuanced discussion of the approach stating how it supports scientific discovery, what is novel relative to the earlier studies, and what are the open challenges. Overall the approach is technically solid and well motivated, and eLife readers will benefit from such a discussion in the context of the new and extensive data set.
Regarding biological significance, the reviewers all agree that the authors do not clearly address how their results advance our understanding of the biology, nor provide any novel biological insights. It is very important that any conceptual novelty is highlighted in the revised manuscript.
Other serious concerns include limitations of the approach with respect to predictions (reviewer 2, please address this in the response and in the manuscript), and lack of a convincing statistical validation (e.g., frequentist coverage of the predictions in a realistic simulation study – where data are generated using parameters estimated from the real data). We would like to see these addressed in the revision as well.
Overall, this is a solid piece of work, and we look forward to receiving a revision.
Reviewer #1:
The authors applied belief propagation to data from systematic pharmacological perturbations and generated network models of signaling in melanoma cells. The authors used these models to predict cellular responses to untested drug perturbations. Simulating the models, the authors made predictions for effective combinations of perturbation. The authors argue that one of these predictions, that cotargeting cMyc with MEK or RAF is synergistic, is non trivial and go on to demonstrate experimentally a synergistic effect of cotreating.
Conceptual and highlevel feedback:
In my opinion, the most significant and general result of this work is showing the promise of genomicallyinformed preclinical trials and providing a concrete example of a synergistic pharmacological treatment inhibiting the cell division of melanoma cells. The experimental design and the network inference methodology seem careful and well executed. The results appear convincing. On the downside, I do not see much conceptual novelty or results contributing to our understanding of the biology. It is also not clear to me that the initial synergistic effect of cotargeting, which is interesting and explored in details, will actually prevent the development of drugresistance over the long term in a clinical application. I invite the authors to emphasize any conceptual advances in their work and new biological conclusions that I may be missing. I think that conceptual advances would contribute the most to elevating the significance and general interest in this work.
[Minor comments not shown.]
Reviewer #2:
In this study, the authors use a belief propagation (BP) algorithm to construct a network model of nearly 100 (phospho)proteins and several phenotypes in response to single and paired drug treatments. The network models use a differential equation framework, so they can be simulated to predict the response of the observed variables (proteins, phenotypes, etc.) to new perturbations and arbitrary combinations of perturbations. Simulations nominate cMyc, which was not perturbed in the original drug screen, due to its in silico effect on G1 arrest. cMyc is targeted indirectly with JQ1, and JQ1 in combination with predicted complementary inhibitors (MEKi and RAFi) induce strong changes in viability and G1 arrest.
The modeling approach is inherently limited to make predictions about proteins on the protein array in the perturbation screen and edges represent indirect effects, but the experimental/computational technique – in particular the focus on quantitative, predictive instead of descriptive networks – is powerful when there is sufficient coverage on the protein array. The algorithmic innovations have already been published (BP algorithm and preliminary SkMel133 analysis in Molinelli 2013; including prior knowledge in Miller 2013). This manuscript expands the (phospho)protein antibodies (16 to 138), perturbation conditions (44 to 89), and phenotypes (1 to 5) modeled in RAFiresistant SkMel133 cells and evaluates new predicted pairs of inhibitors. It also introduces the preprocessing PERA step, which improves the inclusion of prior pathway knowledge by making it more systematic. The prior implementation – using soft constraints and rewarding edges with prior support instead of penalizing those without – is wellmotivated given that interactions in pathway databases are derived from different cell types.
1) The criteria used to select cMyc and G1 arrest instead of other phenotypes and topranked in silico perturbations are not made explicit, making it difficult to assess how well the method could predict additional novel perturbations. Figure 5–figure supplement 2 shows that the top perturbation for all phenotypes is not novel but rather one of the activity nodes for the originally screened drugs. cMyc does cause a substantial change to G1 arrest in silico, but so do other novel perturbations like p38 on viability, STAT3pY705 on S arrest, etc. Furthermore, it appears that there was no preference given to combinations of drugs that are predicted to be synergistic (similar to Miller 2013) or take effect at the lowest possible doses. Indeed, the top viability perturbations – PKCi and CDK4i – only reduced viability at high doses.
2) Is the RAFiresistance phenotype of the SkMel133 cell line recapitulated in the computational model? bRAF or aBRAF are top 10 predicted perturbations for all phenotypes except G2 arrest.
3) Are the observed changes in G1 arrest under inhibitor treatment statistically significant?
4) What was done to verify the identity of the SkMel133 cell line?
5) Edge frequencies are one way to assess BP stability, and Figure 4–figure supplement 2 shows that there are few stable edges. How much do the phenotype predictions vary across the 4000 simulations?
[Minor comments not shown.]
https://doi.org/10.7554/eLife.04640.030Author response
The editorial team has identified two areas of major concern: (1) novelty of the approach and (2) biological significance.
We included two new sections to address these concerns (please see below).
Regarding novelty of the approach, a portion of the most interesting methodological innovations has been already published by a subset of the authors. The preprocessing PERA step is new, and so are a number of other technical advances. In this paper, the authors extended the scope of the analysis substantially (phosphoprotein antibodies from 16 to 138, perturbation conditions from 44 to 89, and phenotypes from 1 to 5), and we appreciate that this extension is critical for biologists to adopt such approaches. However, the authors fail to clearly present the innovations over previous work. The authors should make an effort, perhaps in a dedicated subsection under Methods, to offer a nuanced discussion of the approach stating how it supports scientific discovery, what is novel relative to the earlier studies, and what are the open challenges. Overall the approach is technically solid and wellmotivated, and eLife readers will benefit from such a discussion in the context of the new and extensive data set.
We included a detailed discussion of the methodological innovations, open challenges and significance of the method for scientific discovery. This is a separate section titled “Advances in the perturbation biology method” in Methods.
Regarding biological significance, the reviewers all agree that the authors do not clearly address how their results advance our understanding of the biology, nor provide any novel biological insights. It is very important for eLife that any conceptual novelty is highlighted in the revised manuscript.
We rewrote a significant portion of the Discussion to better address the significance of our biological findings and their importance for cancer research, melanoma and drug resistance. In this section, we summarized:
A) The previous attempts (drug combinations) to overcome the drug resistance in melanoma.
B) Why targeting cMyc and the ERK pathway is potentially more beneficial than other combinatorial interventions in melanoma.
C) A more generalized view of overcoming drug resistance by cotargeting upstream genomic aberrations and downstream signaling proteins, which are connected to diverse pathways.
D) The details of our findings in terms of synergistic response to BET bromodomain inhibitors and ERK pathway inhibitors.
E) Potential clinical applications of BET bromodomain inhibitors either as a single agent or in combination with kinase inhibitors (e.g., RAFinhibitor) in cancer therapy.
Other serious concerns include limitations of the approach with respect to predictions (reviewer 2, please address this in the response and in the manuscript), and lack of a convincing statistical validation (e.g., frequentist coverage of the predictions in a realistic simulation study – where data are generated using parameters estimated from the real data). We would like to see these addressed in the revision as well.
For limitations of the approach with respect to predictions, please see the response to reviewer 2.
For statistical validation, we computed the frequentist coverage probabilities of our predictions (Figure 3, Figure 3–figure supplement 1 and Supplementary file 1E) for all cross validation predictions (CVs). First, we expanded the leave11out crossvalidation (CV) simulations from 2 (AKTi and RAFi in Figure 3) to all possible conditions (2+10 calculations, Figure 3–figure supplement 1). Next, for each crossvalidation/simulation set, we provided the frequentist coverage probabilities in 95%, 90% and 70% confidence intervals (CI) (Supplementary file 1F). The statistical validations show that in 10 out 12 of leave11out CVs, the frequentist coverage is > 95%. For CV that involves the prediction of response to HDACi, the frequentist coverage probability is 93% for the 95% CI. We concluded that this is also a highquality prediction given the level of ambition in our cross validation scheme. However, the frequentist coverage for CDK4i predictions is below 95% (88%). Indeed, the correlation coefficient between the simulated values and the hidden experimental conditions was also ∼0.70 for CDK4i, while the correlation coefficients are as high as 0.88 for other inhibitors. This is most probably due to the fact that the CDK4 inhibitor was highly toxic compared to the other compounds. The number of viable cells under CDK4i was lower than those in other perturbation conditions. Therefore, the prediction quality for CDK4i crossvalidation was lower compared to the other CVs. Although we state that a predictive power of 0.70 is still sufficient to make useful predictions, this demonstrates the importance of a high quality training dataset for reaching excellent predictive power such as in all other conditions.
Overall, this is a solid piece of work, and we look forward to receiving a revision.
Reviewer #1:
Conceptual and highlevel feedback:
In my opinion, the most significant and general result of this work is showing the promise of genomicallyinformed preclinical trials and providing a concrete example of a synergistic pharmacological treatment inhibiting the cell division of melanoma cells. The experimental design and the network inference methodology seem careful and well executed. The results appear convincing. On the downside, I do not see much conceptual novelty or results contributing to our understanding of the biology. It is also not clear to me that the initial synergistic effect of cotargeting, which is interesting and explored in details, will actually prevent the development of drugresistance over the long term in a clinical application. I invite the authors to emphasize any conceptual advances in their work and new biological conclusions that I may be missing. I think that conceptual advances would contribute the most to elevating the significance and general interest in this work.
We thank the reviewer for this suggestion and we believe this will certainly increase the significance of our study. As stated above, we have included two separate sections to address this issue: (1) we updated the Discussion section and described in detail the biological significance of our findings and the importance of the discovered combination strategy for RAFi resistance problem in melanoma biology; (2) we included another section to the Methods, where we discuss in detail the methodological improvements in this paper, how it can accelerate the biological discoveries and what the current open challenges are.
Reviewer #2:
1) The criteria used to select cMyc and G1 arrest instead of other phenotypes and topranked in silico perturbations are not made explicit, making it difficult to assess how well the method could predict additional novel perturbations. Figure 5–figure supplement 2 shows that the top perturbation for all phenotypes is not novel but rather one of the activity nodes for the originally screened drugs. cMyc does cause a substantial change to G1 arrest in silico, but so do other novel perturbations like p38 on viability, STAT3pY705 on S arrest, etc. Furthermore, it appears that there was no preference given to combinations of drugs that are predicted to be synergistic (similar to Miller 2013) or take effect at the lowest possible doses. Indeed, the top viability perturbations – PKCi and CDK4i – only reduced viability at high doses.
We ranked the perturbation conditions based on their predicted influence on each phenotype in the model. We selected the top ranked perturbations as candidates for further experimental validation. One of the strengths of our method is that there is not a single but a set of useful (top ranked) predicted combinatorial perturbations in a given context (e.g., RAFi resistant melanoma cells). Each of these predictions can be tested using targeted drugs or other systems such as shRNA, CRISPRCas. We believe that these predictions will be a very useful resource for experimental biology groups that aim to test drug combinations in melanoma. We now have provided the predicted response to all the combinatorial perturbations on our webpage. Biologists can explore the predicted responses to > 70000 combinatorial perturbations (combinations of different doses, strengths etc.) for five phenotypes. We also provide all of the model solutions so that they can predict response to alternative combination schemes (e.g., high order combinations, time staggered combinations etc.) using any type of program that can execute ODEs.
In this study, we chose to experimentally test one of the most promising predicted combinations: cMyc/ERKpathway (MEKi or RAFi) and G1 arrest. Reasons for selecting this combination include: (1) we predicted a substantial synergy for these combinations and these are the top responding targets for G1arrest; (2) the combinations involve targets (e.g., ERK pathway) that are known targets for treatment of melanoma (3) induction of G1arrest is a promising strategy for controlling melanoma growth (Solit et al., Nature, 439, 358362); and various genes and cell cycle pathways associated with G1arrest are involved in melanoma progression and RAFinhibitor resistance (Rother and Jones, 2009; Xing et al., 2012). Such genes involve RB, CDKN2A & CDKN2B (deleted in SkMel133 cells), CDK4 and many others (4) the predicted targets are potentially druggable, albeit indirectly in the case of cMyc. We provide further details related to other predicted combinations below.
There are certainly other very promising predictions, such as targeting pcJUN (Using JNK inhibitors) and ERK pathway for G1arrest.
Some of the predicted perturbations/targets (activity nodes) were already included in the experimental training set. Combinatorial perturbations that involve such trivial and other novel/nontrivial targets should still be tested in the laboratory. One example is the HDAC inhibitor that induces a strong G2arrest. It will be interesting to search for useful combinations involving HDACi. Some predicted targets are not “activity nodes” (e.g., STAT3pY705) and are less interesting than others since they are highly related to the activity nodes (i.e., perturbations). For example, STAT3i is already in the experimental set and it reduces both STAT3 phosphorylation (STAT3pY705) and cell viability. Therefore it was not an immediate choice for us. Interestingly p53 and aMDM2 have very similar G2arrest responses (Figure 5–figure supplement 2). Note that MDM2 is a p53specific ubiquitin ligase and controls p53 stability.
Some of the predicted targets are not directly druggable with targeted compounds (e.g., most total level targets, such as p38/MAPK14, p27/Kip1, CRAF, TSC2, caveolin). New gene editing technologies such as CRISPR may be useful to experimentally perturb these nodes and test our predictions.
We provide a very brief summary of the above in the revised manuscript (Results, last paragraph). We also included a table that lists the top predictions (Supplementary file 1).
2) Is the RAFiresistance phenotype of the SkMel133 cell line recapitulated in the computational model? bRAF or aBRAF are top 10 predicted perturbations for all phenotypes except G2 arrest.
The RAFinhibitor resistance was recapitulated in the computational models. A good measure of the RAFinhibitor resistance in the models is the response to perturbation on the aBRAF but not the BRAF. The BRAF node corresponds to the total protein level of the BRAF protein and aBRAF is an activity node that couples the effect of the RAF inhibitor to rest of the network.
aBRAF is in the top 10 predictions only for G1arrest. The experimental phenotypic response data (Figure 2C) show a slight G1arrest response to MEK and RAF inhibitors (Note that RAFi/MEKi combination does not lead to a drastic increase G1arrest. Also see the strong G2arrest response to HDACi, which was recapitulated in the models; Figure 5–figure supplement 2). These effects are recapitulated in the models as aBRAF is in top 10 for the G1arrest prediction. However, as can be seen in Figure 6A, the response to combined targeting of cMyc, CyclinD1, aMEK, aBRAF and pcJun is higher than the single agent perturbations. aBRAF was not in the top10 predictions for other phenotypes.
The total BRAF level is in the top 10 for cell viability, Sarrest and G2M nodes but the predictions do not suggest that BRAF has a drastic effect on any of the phenotypes (compared to the top few predictions). According to the vast literature on RAF inhibitors and melanoma, the total BRAF level usually does not immediately change in response to RAF inhibition and is not a direct measure of RAFi resistance, as RAFi inhibits phosphorylation of MEK by mutated BRAF protein. However, it is very interesting that BRAF (and CRAF) levels are related to cell viability and some of the cell cycle phenotypes. It is reported that the longterm changes in expression levels and heterodimerization of BRAF and CRAF have nontrivial implications in the emergence of resistance to RAF inhibitors (See Figure 2 in Lito, Rosen and Solit, 2013, Tumor adaptation and resistance to RAF inhibitors, Nature Medicine). Due to the transactivation of the wildtype BRAF and CRAF and RAF isoform switching in response to RAF inhibition, the relation between the (B/C)RAF concentration and drug resistance is not linear. Experimentally testing the effect of changes in BRAF/CRAF on RAF inhibitorinduced phenotypes and RAFi sensitivity will be very interesting. On the other hand, we did not perform experimental tests related to predictions on RAF levels. In the current manuscript, it would not be more than a speculation to relate (B/C)RAF levels and drug resistance.
3) Are the observed changes in G1 arrest under inhibitor treatment statistically significant?
Cell cycle arrest experiments (flow cytometry) were run in triplicates. Although we refrained from running a ttest with three data points, we provided standard error bars in Figure 6D. The error bars indicate that the G1 arrest response for the combination is significantly stronger than those for other conditions (no drug and single drug perturbations).
4) What was done to verify the identity of the SkMel133 cell line?
Our collaborators and coauthors, David Solit and Christine Pratilas, verified the identity of the SkMel133 and a panel of other melanoma cell lines. They checked for mislabeling, contamination and misidentification using a multiplexed PCR/mass spectrometry (MS)–based genetic fingerprinting assay run on the HapMap Sequenom platform. The method is described in: Janakiraman et al., Genomic and Biological Characterization of Exon 4 KRAS Mutations in Human Cancer, Cancer Res. 2010 Jul 15;70(14):590111. We now state this in the Methods section of the revised manuscript.
5) Edge frequencies are one way to assess BP stability, and Figure 4–figure supplement 2 shows that there are few stable edges. How much do the phenotype predictions vary across the 4000 simulations?
We now provide the distributions of all the predicted responses over 4000 solutions on our webpage (http://130.211.159.201/pertbio/). The statistical information for each distribution is also provided on the webpage.
https://doi.org/10.7554/eLife.04640.031Article and author information
Author details
Funding
National Institutes of Health (NIH) (The National Resource for Network Biology grant, GM 103504)
 Anil Korkut
 Emek Demir
 Özgün Babur
 Chris Sander
National Human Genome Research Institute (NHGRI) (The Research Resource for Biological Pathways grant, U41 HG006623)
 Emek Demir
 Özgün Babur
 Chris Sander
Melanoma Research Alliance (MRA) (Established Investigator Award)
 Anil Korkut
 Weiqing Wang
 Chris Sander
National Institutes of Health (NIH) (Center for Cancer Systems Biology grant, U54 CA148967)
 Anil Korkut
 Weiqing Wang
 Chris Sander
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank NP Gauthier, P Kaushik, ML Miller, DS Marks, A Braunstein, A Pagnani, R Zecchina, M Cokol, and N Rosen for helpful discussions. We thank G Mills for critical reading of the manuscript, helpful discussions, and help with the acquisition of the RPPA data. This work was funded in part by the Center for Cancer Systems Biology grant U54 CA148967 (NIH), the National Resource for Network Biology grant GM103504 (NIH), the Research Resource for Biological Pathways grant U41 HG006623 (NHGRI), and a Melanoma Research Alliance Established Investigator Award.
Reviewing Editor
 Edoardo M Airoldi, Harvard University, United States
Publication history
 Received: September 10, 2014
 Accepted: July 7, 2015
 Version of Record published: August 18, 2015 (version 1)
Copyright
© 2015, Korkut et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 4,421
 Page views

 877
 Downloads

 45
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.