article: We present a powerful experimental-computational technology for inferring network models that predict the response of cells to perturbations, and that may be useful in the design of combinatorial therapy against cancer. The experiments are systematic series of perturbations of cancer cell lines by targeted drugs, singly or in combination. The response to perturbation is quantified in terms of relative changes in the measured levels of proteins, phospho-proteins and cellular phenotypes such as viability. Computational network models are derived de novo, i. e. , without prior knowledge of signaling pathways, and are based on simple non-linear differential equations. The prohibitively large solution space of all possible network models is explored efficiently using a probabilistic algorithm, Belief Propagation (BP), which is three orders of magnitude faster than standard Monte Carlo methods. Explicit executable models are derived for a set of perturbation experiments in SKMEL-133 melanoma cell lines, which are resistant to the therapeutically important inhibitor of RAF kinase. The resulting network models reproduce and extend known pathway biology. They empower potential discoveries of new molecular interactions and predict efficacious novel drug perturbations, such as the inhibition of PLK1, which is verified experimentally. This technology is suitable for application to larger systems in diverse areas of molecular biology.
  Abnormal biomolecular information flow as a result of genetic or epigenetic alterations may lead to tumorigenic transformation and malignancy and is classically modeled as changes in signaling pathways [1]. Targeted anti-cancer drugs, which bind and inhibit specific components of aberrant signaling pathways, are a promising alternative to conventional chemotherapy, with recent successes in melanoma (RAF inhibitor) [2] and prostate cancer (AR inhibitor) [3], [4] following in the footsteps of the pioneering BCR-ABL inhibitor Imatinib [5] and EGFR inhibitors Gefitinib and Erlotinib [6], [7], [8]. Combinations of targeted anticancer drugs hold considerable promise because of the emergence of resistance to initially successful single agents and the highly robust nature of the signaling pathways with multiple feedback mechanisms [9]. High throughput measurements on response profiles of living cells to multiple perturbations such as drug combinations provide a rich set of information to construct quantitative cell biology models. In this paper, we construct context specific de novo mathematical models of signaling pathways through the use of systematic paired perturbation experiments and network inference algorithms. Such network models provide insight into mechanistic details of signaling pathways, predict the response of cellular systems to multiple perturbations beyond those from which models are derived, and guide the design of perturbations for a desired response. Previous mathematical models of molecular signaling in cells have been effective in modeling pathways and enhancing drug discovery [10], [11], [12], [13], [14], [15], [16], [17], [18], [19]. Techniques for network modeling of signaling pathways span a wide spectrum of complexity. Detailed chemical kinetics and spatiotemporal models [17], [20], [21] can provide mechanistic explanations of observed behavior, but are often incompletely parameterized, needing tens or hundreds of sensitive parameters for medium-sized systems. Moreover, such models may not be valid in biological contexts that differ substantially from dilute solution chemistry. On the other end of the spectrum, pattern matching or machine learning models such as neural networks and correlation-based models such as maximum entropy [22] can accurately provide purely data-driven models of signaling. However, such methods have limited power to explain mechanistic details and, in most cases, are insufficient for quantitative predictions of system behavior in conditions beyond those from which the models are derived. We take a unique modeling approach to construct context specific, de novo and predictive network models of signaling pathways from drug perturbation data (Figure 1). Here, de novo means that network inference is done without depending on known molecular interactions extracted from literature or pathway databases, which do not account for biological context. This approach also emphasizes context specificity since it relies on rich experimental data from a single biological context as its training set. The models are constructed through parameterization of a simple model equation, which has been used in other network modeling approaches [13], [23], [24], [25], [26], [27]. The model equations contain parameters that are mechanistically descriptive of direct or indirect interactions in the system. Finally, the models are computationally predictive of cell-type specific response to new drug perturbations and their combinations. We expect that this conceptual framework and the technical advances in network inference will empower the community to identify unique drug targets and combinations that are particularly efficacious within specific disease contexts. De novo construction of signaling network models at scales relevant to problems related to complex biological phenomena such as cancer has long been a challenge in system biology. Thus, quantitative models of protein signaling pathways are typically constructed on the basis of existing prior knowledge from literature searches [16], [17] and interaction databases [19], [28], [29]. However, different cancer contexts have unique genetic and proteomic alterations to normal protein signaling. For example, distinct mutations in effector proteins of the PI3K pathway are oncogenic in unique ways, i. e. , they lead to distinctly context dependent functional consequences in different cancer types, subtypes and patients [30]. The method we introduce here is capable of inferring parameterized network models with prior knowledge, yet the results in this study do not include prior knowledge. An advantage of de novo inference is independence from prior knowledge interactions that may be incorrect or incomplete in a particular biological context. The largest obstacle to de novo model construction is the combinatorial explosion in the number of possible network models, which defines the solution space [31]. Model inference problems of this type are NP-hard [32]. The number of possible configurations for a model with N nodes, and K possible values for each parameter grows super-exponentially as. We have previously described a method named CoPIA for de novo construction of dynamic nonlinear network models from perturbation data [33]. CoPIA is based on the combined use of a Monte Carlo stochastic search algorithm, which is used to search the network configurations and an efficient gradient descent algorithm [34] for quantitative parameter optimization. However, without algorithmic improvements, such Monte Carlo based methods are limited to modeling fewer than approximately 15–20 biological entities [33]. Increasing the scale of network (or pathway) models of cellular signaling processes to levels sufficient to describe complex biological problems in quantitative detail is therefore extremely challenging and has been approached with a diversity of methods [35]. An ingenious, two-step approach to deal with network inference in larger systems is based on first calculating probability distributions for each possible interaction in the model and then computing distinct solutions by sampling these probability distributions. For this purpose, we employ a probability model of network configurations inspired from statistical physics principles. Following a set of approximations to simplify the probability model, we apply a custom adaptation of an iterative algorithm called Belief Propagation (BP). BP involves local optimization updates to probability distributions of individual model parameters that converge to a stable set of probability distributions, which collectively describe a set of good network model solutions [36], [37], [38]. BP has been applied to various complex inference problems, some of them NP-Hard such as K-SAT [38] and graph coloring [39]. BP has garnered some attention in biological network inference [40], [41] and parameter estimation [42], [43]. Here, we tailor the BP algorithm to large-scale perturbation data that is capable of increasing the scope of the models to hundreds of nodes. The result of BP is a set of probability distributions for each model parameter, often referred to as marginal probability distributions, or ‘marginals’. Each marginal describes the inferred distribution of a particular parameter across a range of high probability solutions. Individual models are created via sampling from these marginals [44]. Consequently, the time-complexity of the problem is strongly reduced, the prohibitive cost from combinatorial complexity is circumvented and, although the method provides only an approximate solution, one obtains useful, non-trivial results. Our algorithmic network pharmacology approach involves four major steps: (i) perturbation experiments with combinations of targeted compounds; (ii) high-throughput quantitative measurements of proteomic changes (e. g. , reverse phase protein arrays or mass spectrometry) and phenotypic changes (e. g. , cell viability or apoptosis); (iii) inference of quantitative network models of protein signaling that explain and link these changes; and (iv) use of the network models to predict cellular and molecular responses to diverse perturbations, beyond the conditions on which the network models are derived. In this work, we adapt BP to construct quantitative network models of signaling pathways from systematic perturbation experiments. We evaluate the speed and accuracy of BP on toy data generated from biologically inspired network structures. The inference on this toy data reveals that BP offers a significant improvement in computational efficiency compared to traditional Monte Carlo simulations without a sacrifice in accuracy. Furthermore, we construct network models of signaling in a RAF inhibitor resistant melanoma cell line (SKMEL-133), which has the BRAFV600E mutation [45], [46]. The models are predictive of both the proteomic and phenotypic response to drug combinations. Model simulations successfully predict the phenotypic response profiles of SKMEL-133 cells to novel drug targets. With the introduction of many novel targeted drugs and patient specific genomic profiling, this network pharmacology approach aims to provide an effective tool to develop individualized combination therapies against multiple cancer forms.
  The success of BP depends on whether or not the simulations of the models taken from BP are quantitatively predictive of cellular response to drug combinations. It is also useful, as an exercise, to evaluate the overall performance of the BP algorithm on data sets engineered from completely known networks. With such toy datasets we achieve the following: (i) demonstrate that BP converges quickly and correctly; (ii) compare BP network models to a known data-generating network; and (iii) evaluate performance in biologically realistic conditions of noisy data from sparse perturbations. The synthetic data is generated without the assumptions used in the formulation of BP, and therefore serves as a reasonable test of the sensitivity of the BP method to those assumptions. See methods for more information on how the toy data is generated. The probabilistic nature of the BP algorithm is the key feature that enables de novo inference on large and complex problems in cell biology, e. g. , signaling processes involving more than a hundred molecular and phenotypic variables, which previous methods could not reach. Given this opportunity, we apply our network pharmacology approach to SKMEL-133, a melanoma cell line resistant to the RAF inhibitor (Vemurafenib, PLX4032), which inhibits the BRAFV600E mutant protein kinase more strongly than other RAF proteins.
  We describe a combination of experimental and computational methods, in the field of network pharmacology, to construct quantitative and predictive network models of signaling pathways. The particular contribution is a set of algorithmic advances, which we adapted from statistical physics to infer network models in sizes and complexities not reachable by classical gene-by-gene molecular biology. The necessity of inference of complex network models stems from the fact that classical methods, in which a small number of perturbation experiments lead to functional description of carefully selected sets of genes and gene products is reaching technical limits. High-throughput proteomic and genomic profiling technologies provide much richer and more complex information about cellular responses than can be analyzed by a scientist' s thought processes. At these levels of completeness and detail, predicting changes in physiological attributes from molecular data requires computational modeling and quantitative analysis. Our quantitative network models not only capture already known biological interactions but also nominate novel interactions and may detect complex regulatory mechanisms, such as feedback loops, in specific biological contexts. The quantitative analysis of molecular and cellular behavior in these models provides detailed understanding of the coupling between signaling processes and global cellular behavior. Such understanding is hard to achieve by reductionist approaches that focus on the relation between single or few molecules and cellular processes. Furthermore, we provide a systems biology platform to predict the cellular (i. e. , proteomic and phenotypic) response profiles to multiple perturbations such as those in combinatorial cancer therapies. This investigation constitutes a proof-of-principle of a particular technology: a combination of a network inference algorithm and a technology for perturbing cells and measuring their response. The overall context is network pharmacology, by which we mean the science of using network models to derive and then test effective therapeutic targets and combinations. The particular challenge in cancer biology is the complexity and individual variation of genetic and epigenetic alterations that are plausibly cancer causing, and thus modulate the response to therapy; and, the emergence of resistance to successful targeted therapies, such as EGFR inhibitors in lung cancer or RAF inhibitors in melanoma. In our view, the fairly fragmented gene-by-gene classical methods of molecular biology, while extremely powerful as a reductionist method, are reaching a clear limit. Those methods struggle with effects such as ‘cross-pathway coupling’, ‘multigenic diseases’ or individual variation of response to therapy. More comprehensively quantitative and computationally predictive methods (‘systems biology’) are likely not only to increase our predictive abilities but also save substantial overall effort by computationally testing large numbers of cellular states on a large variety of genetic backgrounds and in exhaustively explored perturbation conditions. The implementation of our BP algorithm enables us to infer much larger network models than are reachable with standard Monte Carlo search methods. First, we use biologically realistic toy data to illustrate the dramatic speed improvement of the BP method over standard MC methods. BP convergence times scale favorably with the size of the models at no measureable loss of accuracy. This is a desirable property for constructing large network models de novo. Even in the biologically realistic conditions of noisy data from sparse perturbation, BP infers most of the dominant interactions in the data-generating network. Although a fraction of the likely BP interactions are false positives, these false positives are consistent with the data and connect structurally correlated nodes. The assumptions inherent in this BP implementation are potential limitations to the overall efficacy and accuracy of our approach. However, assumption 1, while important, does not dramatically affect the results beyond a critical level of discretization. Furthermore, since BP performs as well or better than MC, which employs neither assumption 3 nor 4, we suspect by process of elimination that assumption 2 is the largest contributor to the observed limitations of our current method. For real biological data, the quality of BP' s performance depends on the applicability of the model equations and both the richness and quality of the data used for training. We previously published an MCF7 breast cancer cell line dataset, for which network models were derived with a nested MC search algorithm (CoPIA) [97]. Reassuringly, many of the strongest interactions from CoPIA are recovered by the BP method (Supplementary Text-S4, Figure S10). We have also modeled the same SKMEL-133 dataset with Gaussian Graphical Models, a popular probabilistic model that goes beyond pairwise correlations to distinguish between direct and indirect correlations. A comparison to the BP derived interactions suggests non-trivial overlap in the strongest couplings. A more detailed analysis is available in the supplement (Supplementary Text –S5, Figures S11 and S12). Finally, we have also compared BP against both prior knowledge models and Bayesian network models inferred with Bayesian inference [98]. That analysis indicates superior performance for the BP derived models for at least this kind of data and this style of model (See Supplemental Text S7 and Figure S17). We are able to capture the known interactions in MAPK and AKT/PI3K canonical pathways through de novo network inference. Moreover, we predict potential interactions, either direct or indirect, such as a bidirectional interaction between Rb and MAPK pathways and a potential feedback loop from PI3K and MAPK pathways acting on SRC activity through inhibitory phosphorylation at SRCpY527. We also quantitatively predict the effect on cell viability from inhibition of novel targets via explicit in silico simulation. The in silico perturbations of Cyclin B1 and PLK1 lead to comparable reduction in cell viability. As neither Cyclin B1/PLK1 nor their direct regulators are inhibited for the model training dataset, this result is a genuine and nontrivial prediction. We test and validate the PLK1 prediction by measuring the cell viability after treating SKMEL-133 cells with a potent and selective PLK1 inhibitor. Indeed, the PLK1 inhibitor reduces cellular growth significantly in SKMEL-133 cells even when treated with the drug at nanomolar concentrations (Figure 9). Moreover, none of the in silico perturbations of the MAPK pathway lead to a significant change in cell viability, which is consistent with the experimental responses to perturbations acting on MEK and BRAFV600E and RAF inhibitor resistance of SKMEL-133 cells. We are aware that drugs do not usually have a single specific target. While this is potentially problematic for modeling and simulating drug effects, we find that the inference is driven largely by the correlations in the data arising from the effect of perturbations on the overall system. We estimate that BP is most sensitive to strong and untrue assumptions about the direct effects of a given drug. The use of activity nodes is an indirect way of dealing with off-target drug effects. These nodes, which are perturbed but not measured, represent the coupling of a drug to the rest of the system. Such coupling includes both specific and off-target effects. BP may infer interactions from these activity nodes to any number other measured nodes thus simultaneously inferring targets for the particular drug used. Drug specificity also affects the predictive power of the resulting models. All effects from the AKT inhibitor, for example, are assigned to single inhibition of the aAKT node, even if other targets of the AKT inhibitor are partially responsible for the measured outcomes. Therefore, any simulation-based prediction regarding perturbation of the aAKT node will be tethered to the off-target effects of the AKT inhibitor used in training. The complete solution of the drug specificity problem requires a more comprehensive and systematic analysis to determine the effect of off-target effects on quality and predictive power of inferred models. The leave-k-out analysis reveals few outlier points, where the predicted response profiles largely disagrees with test data. These (mis) predictions fall into two major categories. In the first category, (mis) predictions arise due to measurements with very low signal to noise ratios and high experimental uncertainties. The (mis) predicted S6pS240 levels and cell viability phenotypes in some of the perturbation conditions fall into this category. Note that the models are trained and simulated after logarithmic conversion (i. e. , xi = log2-ratios of measured signals), which exaggerates the errors for low signals. In linear space, no such outliers are observed suggesting that those outliers can be considered as an artifact of our analysis in logarithmic space. In the second category, (mis) predictions arise due to insufficient experimental constraints. (Mis) prediction of AKTpS473 levels after mTOR inhibition falls into this category. mTOR inhibition leads to an increase in AKT phosphorylation possibly due to the disruption of a feedback loop. In our current analysis, mTOR inhibition and steady state measurements on AKT phosphorylation are the sole experimental inputs to detect the changes in this feedback mechanism. When the mTOR inhibition is withheld in leave-k-out tests, experimental constraints become insufficient to describe the regulation of AKT phosphorylation, thus leading to the (mis) predictions. A systematic experimental design to enrich the perturbations, better characterization of the AKT phosphorylation dynamic range and richer proteomic measurements on this part of signaling pathways are possible ways to improve the quality of these particular predictions. In general, careful optimization of perturbation conditions (drugs and combinations) and observations (protein arrays, mass spectrometry target list), within available resources, would significantly enhance the predictive power of this approach. Predicting the effect of drug combinations is highly challenging and has been the subject of many studies, e. g. , [99], [100]. When a large repertoire of targeted drugs are available for screening, the search space for useful combinations of two or more drugs is combinatorially complex. Moreover, therapeutically promising drug combinations are not limited to simultaneously introduced perturbations. Potentially useful drug combinations may consist of combinations of relative doses, for two, three and four drugs applied simultaneously or sequentially after well-defined time intervals. Identification of such complicated combinations through experimental screening tests is prohibitively cumbersome. Here, we provide a potential solution to this problem in that in silico screens using predictive network models can cover a large space of possible drug combinations. The predictive power of the models derived here is apparent from the reasonable accuracy of predicting the results of withheld experiments using a leave-k-out cross validation. Network models inferred from perturbation experiments can thus be used to predict the responses to novel combinations. One can enlarge the search space for drug treatments from a few hundred experimentally screened combinations to tens of thousands of computationally tested combinations and guide subsequent, highly efficient experimental screens of the top predicted candidates. Beyond the computational power of well-constrained and robustly derived network models, one may expect to achieve a conceptual understanding of the principles of epistasis of drug effects and the mechanisms of resistance to targeted therapeutics. For example, the initial system response to drug intervention on the scale of minutes to days may be indicative of subsequent epigenetic and genetic changes in a population of treated cells that represent the long-term and hard-to-treat emergence of drug resistance. In this context, reliable dynamic network modeling may be an excellent guide to strategies for blocking the emergence of resistance in the first place. The pre-clinical consequence is the selection of combinations of therapeutic interventions that not only are effective in slowing the proliferation or promoting the elimination of cancer cells, but also counteract resistance to otherwise effective treatments, such as RAF inhibitors in melanoma or AR inhibitors in advanced prostate cancer, in clinical trials. While the approximate solutions to the problem of network inference deliver interpretable biological results, there is much room for improving the power and information value of the method. For example, when time-dependent response measurements after perturbation are available, one can use analogous, extended algorithms to infer probability distributions over (time-independent) interaction parameters W that describe the behavior of a time-dependent system. Another extension of the formalism is the use of more complicated forms of the differential equations, such as those of enzyme kinetics. Also, even in the current approximation, careful design of experiments selecting a minimal set of maximally informative perturbation conditions would increase the efficiency of this experimental-theoretical approach. A straightforward and powerful extension is the systematic use of prior information in the form of directed interactions adapted from the current scientific literature or pathway databases. Such prior information is easily incorporated as a set of additional constraints on the probability distributions of W in Equation 3. This has been utilized in our related work on modeling dedifferentiated liposarcoma drug combinations [43]. We are also actively pursuing a method for the systematic inclusion of prior knowledge interactions from curated databases. On the experimental side, measurement of richer phenotypic attributes of cells, such as apoptosis or cell cycle arrest, as well as markers of differentiation states would greatly increase the predictive power by providing more links between molecular and phenotypic quantities. The network pharmacology approach described here provides a strong tool for a system level description of signaling events in cancer cells. Moreover, it presents a step forward in quantitative prediction of responses of cancer cells to drug perturbations. Beyond cancer biology, there is no reason to believe that the proposed technology cannot be used to derive accurate quantitative and predictive network models for biological cellular systems in general, provided sufficiently diverse experimental perturbations and sufficiently rich readouts are accessible. In this way, we hope to extend the power of classical molecular biology to a broad spectrum of cellular systems with targeted, and possibly clinical, applications.
summary: Drugs that target specific effects of signaling proteins are promising agents for treating cancer. One of the many obstacles facing optimal drug design is inadequate quantitative understanding of the coordinated interactions between signaling proteins. De novo model inference of network or pathway models refers to the algorithmic construction of mathematical predictive models from experimental data without dependence on prior knowledge. De novo inference is difficult because of the prohibitively large number of possible sets of interactions that may or may not be consistent with observations. Our new method overcomes this difficulty by adapting a method from statistical physics, called Belief Propagation, which first calculates probabilistically the most likely interactions in the vast space of all possible solutions, then derives a set of individual, highly probable solutions in the form of executable models. In this paper, we test this method on artificial data and then apply it to model signaling pathways in a BRAF-mutant melanoma cancer cell line based on a large set of rich output measurements from a systematic set of perturbation experiments using drug combinations. Our results are in agreement with established biological knowledge, predict novel interactions, and predict efficacious drug targets that are specific to the experimental cell line and potentially to related tumors. The method has the potential, with sufficient systematic perturbation data, to model, de novo and quantitatively, the effects of hundreds of proteins on cellular responses, on a scale that is currently unreachable in diverse areas of cell biology. In a disease context, the method is applicable to the computational design of novel combination drug treatments.
title: Perturbation Biology: Inferring Signaling Networks in Cellular Systems