Title: Models from experiments: combinatorial drug perturbations of cancer cells
Abstract: Report2 September 2008Open Access Models from experiments: combinatorial drug perturbations of cancer cells Sven Nelander Corresponding Author Sven Nelander Computational Biology center, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Weiqing Wang Weiqing Wang Computational Biology center, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Björn Nilsson Björn Nilsson Department of Clinical Genetics, Lund University Hospital, Lund, Sweden Search for more papers by this author Qing-Bai She Qing-Bai She Molecular Pharmacology and Chemistry Program, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Christine Pratilas Christine Pratilas Molecular Pharmacology and Chemistry Program, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Neal Rosen Neal Rosen Molecular Pharmacology and Chemistry Program, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Peter Gennemark Peter Gennemark Department of Mathematical Sciences, University of Gothenburg, Gothenburg, Sweden Search for more papers by this author Chris Sander Chris Sander Computational Biology center, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Sven Nelander Corresponding Author Sven Nelander Computational Biology center, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Weiqing Wang Weiqing Wang Computational Biology center, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Björn Nilsson Björn Nilsson Department of Clinical Genetics, Lund University Hospital, Lund, Sweden Search for more papers by this author Qing-Bai She Qing-Bai She Molecular Pharmacology and Chemistry Program, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Christine Pratilas Christine Pratilas Molecular Pharmacology and Chemistry Program, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Neal Rosen Neal Rosen Molecular Pharmacology and Chemistry Program, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Peter Gennemark Peter Gennemark Department of Mathematical Sciences, University of Gothenburg, Gothenburg, Sweden Search for more papers by this author Chris Sander Chris Sander Computational Biology center, Memorial Sloan-Kettering Cancer Center, New York, NY, USA Search for more papers by this author Author Information Sven Nelander 1, Weiqing Wang1, Björn Nilsson2, Qing-Bai She3, Christine Pratilas3, Neal Rosen3, Peter Gennemark4 and Chris Sander1 1Computational Biology center, Memorial Sloan-Kettering Cancer Center, New York, NY, USA 2Department of Clinical Genetics, Lund University Hospital, Lund, Sweden 3Molecular Pharmacology and Chemistry Program, Memorial Sloan-Kettering Cancer Center, New York, NY, USA 4Department of Mathematical Sciences, University of Gothenburg, Gothenburg, Sweden *Corresponding author. E-mail correspondence reaches the principal authors at multiple-perturbation>at<cbio.mskcc.org Molecular Systems Biology (2008)4:216https://doi.org/10.1038/msb.2008.53 PDFDownload PDF of article text and main figures. ToolsAdd to favoritesDownload CitationsTrack CitationsPermissions ShareFacebookTwitterLinked InMendeleyWechatReddit Figures & Info We present a novel method for deriving network models from molecular profiles of perturbed cellular systems. The network models aim to predict quantitative outcomes of combinatorial perturbations, such as drug pair treatments or multiple genetic alterations. Mathematically, we represent the system by a set of nodes, representing molecular concentrations or cellular processes, a perturbation vector and an interaction matrix. After perturbation, the system evolves in time according to differential equations with built-in nonlinearity, similar to Hopfield networks, capable of representing epistasis and saturation effects. For a particular set of experiments, we derive the interaction matrix by minimizing a composite error function, aiming at accuracy of prediction and simplicity of network structure. To evaluate the predictive potential of the method, we performed 21 drug pair treatment experiments in a human breast cancer cell line (MCF7) with observation of phospho-proteins and cell cycle markers. The best derived network model rediscovered known interactions and contained interesting predictions. Possible applications include the discovery of regulatory interactions, the design of targeted combination therapies and the engineering of molecular biological networks. Introduction Our ability to measure increasingly complete and accurate molecular profiles of living cells motivates new quantitative approaches to cell biology. For example, a key aim of systems biology is to relate changes in molecular behavior to phenotypic consequences. To achieve this aim, computational models of cellular processes are extremely useful, if not essential. Computational models can be used for the analysis of experimental data, for the prediction of outcomes of unseen experiments and for planning interventions designed to modify system behavior. We have developed a particular approach to constructing, optimizing and applying computational models of cellular processes, which we call Combinatorial Perturbation-based Interaction Analysis (CoPIA). The key ingredients of the approach are combinatorial intervention, molecular observation at multiple points, model construction in terms of nonlinear differential equations, optimization of model parameters with simplicity constraints and experimental validation. The power of combinatorial perturbation In molecular biology, a targeted perturbation typically inhibits or activates function of biomolecules, e.g. as a result of drug action, small RNA interference, genetic or epigenetic change (Figure 1). In a single experiment, targeted perturbations can be applied either singly or in combination. Combined perturbation by several agents can be much more informative than that by a single agent, as its effects typically reveal downstream epistasis within the system, such as non-additive synergistic or antagonistic interactions. In addition, a large number of independently informative experiments can be performed if in each experiment a different small set of, e.g. two or three, perturbants is chosen from a larger repertoire. Thus, combinatorial perturbations are potentially powerful investigational tools for extracting information about pathways of molecular interactions in cells (such as A inactivates B, or X and Y are in the same pathway) (Avery and Wasserman, 1992; Kaufman et al, 2005; Kelley and Ideker, 2005; Segre et al, 2005; Yeh et al, 2006; Lehár et al, 2007). Combinatorial perturbations can also be powerful application tools when rationally designed to achieve desired effects. For example, combination of targeted drugs is considered a promising strategy to improve treatment efficacy, reduce off-target effects and/or prevent evolution drug resistance (Borisy et al, 2003; Keith et al, 2005; Komarova and Wodarz, 2005; Chou, 2006). Figure 1.Combinatorial perturbation and multiple input–multiple output (MIMO) models. Upper left: intuitive view of perturbations and their points of action. Small inhibitory RNAs alter gene expression; natural protein ligands and small compounds act, e.g., on receptors, transporters or enzymes. Genetic alterations have diverse functional effects. Perturbations can be natural or investigational. Observations (readouts) typically focus on a phenotype of interest, such as growth or differentiation, and on observation points that are both practical and informative, such as transcripts, protein levels or covalent modifications, e.g. phosphorylation. Upper right: MIMO model. All key system variables are represented as real number variables, the combinatorial perturbations ui as well as their targets, internal variables, observation points and phenotypic outputs yi. Inputs (upper layer, squares) affect the different dynamical variables of the system (circles), some of which can be observed (lower layer, squares). The model represents a processing system that relates the input to the output through the interacting dynamical variables. Representation of coupled perturbations (nonlinear effects) is a key requirement of the modeling method. When rate equations are linear (lower left), perturbation effects will combine additively. However, in a well-parameterized system with nonlinear transfer functions (lower right), epistasis effects will arise, and downstream effects depend on pathway organization. Below: uses of a derived MIMO model. When inputs and outputs are known, a model can be used to infer the internal mechanism of the system (Interpretation). When the inputs and the system are known, the model can be used to predict a response (Prediction). When the system and the desired output are known, the model can serve to design optimal modes of control (Control). Download figure Download PowerPoint With recent advances in molecular technologies—e.g., targeted perturbation by small molecules, full-genome libraries of small RNAs, highly specific antibody assays, massive parallelization and imaging techniques—there is intense interest in the investigational power of multiple perturbation experiments in a variety of biological systems. The inherent complexity of such experiments raises significant challenges in data analysis and an acute need for improving modeling approaches capable of capturing effects such as time-dependent responses, feedback effects and nonlinear couplings. Deriving system models from combinatorial perturbation experiments Computer simulation of pre-defined pathways can be used to predict epistasis effects and explore how pathway organization shapes the perturbation response (Omholt et al, 2000; Segre et al, 2005; Lehár et al, 2007). In many situations however, observational data are provided but the pathway is unknown or only partially known. To solve this problem, our computational modeling approach enables users to construct a complete differential equation model for a system from combinatorial perturbation experiments. In the context of this paper, the system of interest is defined by a particular type of cell, its environment, a time interval of observation and a phenotypic change, such as cell death or growth. The system is further characterized by its points of intervention, such as drug targets, and the points of observation, such as the phosphorylation state of proteins involved in signaling processes (Figure 1). To represent such a system mathematically, we choose network models in which nodes represent molecular concentrations or levels of activity and edges reflect the influence of one node on the time derivative of another. The time evolution of the system is modeled by linear differential equations, modified by a nonlinear transfer function to reflect properties of the system that are not explicitly modeled (Figure 1). We present efficient optimization algorithms to find models that achieve maximum agreement between observation and prediction. Our algorithm is based on a combination of a gradient descent method (to set dynamical parameters) and a Monte Carlo process (to explore alternative network connectivities). We make a software implementation of CoPIA available as platform-independent software (http://cbio.mskcc.org/copia). Testing the predictive power of derived system models We perform combinatorial perturbation experiments in an MCF7 breast cancer cell line to test the modeling framework in the steady-state limit. In this test, we demonstrate how observation of the effects of drug pair perturbations can be exploited to deduce a network model of signaling and phenotype control (reverse engineering of pathways). We use observed molecular state and growth phenotype responses to build predictive models and use these to explain the perturbation–phenotype relationship in terms of coupling between proteins in the EGFR/MAPK and PI3K/AKT pathways. Without using known pathway biology, the resulting model reproduces known regulatory couplings and negative feedback regulation downstream of EGFR and PI3K/AKT/mTOR, and makes predictions about possible roles of PKC-δ and eIF4E in the control of MAPK signaling and G1 arrest in MCF cells. We conclude that CoPIA may be of interest as a broadly applicable tool to construct models, discover regulatory interactions and predict cellular responses. For instance, researchers can measure a set of protein phosphorylation responses to drug combinations and use the method to automatically construct network models that predict the response to novel drug combinations. Application of this methodology to time-dependent experimental observations would extend this predictive capability to the regimen of time-dependent, rationally designed combinatorial therapy. Results Modeling the effects of combinatorial perturbations Multiple input–multiple output models State space representation is commonly used in mathematical modeling of input–output behavior in natural systems. In this representation, the time behavior of the system state is described by a first-order differential equation where the vector y(t) represents state variables (the activities of the system's components), the vector u(t) represents perturbations (external influences on the components) and f is a linear or nonlinear transfer function (de Jong, 2002). For example, y(t) can be the abundances of specific mRNAs or proteins, whereas u(t) can be the concentrations of different chemical compounds to which the cells are exposed (Figure 1). In essence, state space models relate a system's input to its output. State space models with multiple inputs–outputs (that is, y and u have more than one coordinate) are called multiple input–multiple output (MIMO) models. Linear MIMO models When f is a linear function of y and u, the above model is called a linear MIMO model. The mathematical properties of linear MIMO models are well known (Ljung, 1986) and such models have been applied to many biological problems, for example, the construction of transcriptional network models (Tegner et al, 2003; Xiong et al, 2004; di Bernardo et al, 2005). Nevertheless, linear models have a limitation in that they can only model uncoupled perturbation effects (linear dose–response relationships), whereas nonlinear effects (coupled perturbation effects) are ignored (Figure 1; 'Model representation'). As a result, linear MIMO models are unable to capture important phenomena that are known to occur in cellular systems, such as saturation effects, switch-like effects and nonlinear interaction phenomena such as genetic epistasis and pharmacological synergism. Nonlinear MIMO models To overcome this limitation, we construct nonlinear MIMO models capable of representing coupled perturbation effects. Previously, other authors have observed that complex gene knockout effects, including epistasis effects, can be predicted in metabolic flux networks where bounds on the reaction rates are introduced (Fell and Small, 1986; Edwards and Palsson, 2000; Segre et al, 2005; Deutscher et al, 2006). Similarly, metabolic systems with Michaelis–Menten kinetics or transcriptional networks with bounds on transcription rates will exhibit epistasis behavior (Omholt et al, 2000; Lehár et al 2007). In the particular case of the MIMO model, we expect more biologically realistic behavior if one replaces the linear transfer function f with a nonlinear transfer function ϕ that imposes bounds on the rates of change of the system. Accordingly, we propose the class of models In this class of models, the matrix wij represents the interactions between the molecules and processes represented by the state variables of the system. (Intuitively, the matrix elements wij can be thought of as a map of the system, in which wij>0 means 'node j activates node i', whereas wij<0 corresponds to inhibition.) Furthermore, αi>0 represents the tendency of the system to return to the initial state (yi=0); βi>0 are constants and ϕi is a transfer function capable of capturing both switch-like behavior and bounded reaction rates. Examples of such functions include sigmoid functions, piece-wise linear approximations of sigmoids or biochemically motivated approximations such as the Hill or Michaelis–Menten equations (Materials and methods). Application of nonlinear MIMO models to combinatorial perturbation experiments We developed computer algorithms to infer nonlinear models of the above type from experimental data, as specified by the best-performing values of the coupling parameters wij and other parameters. As detailed in Materials and methods, the current implementation of our approach consists of the following steps. First, the system of interest is subjected to a set of independent single or multiple target perturbation experiments; and, for each perturbation vector (time-independent instance of u), a readout vector (steady-state instance of y) is recorded. Second, we infer a nonlinear model that best reproduces the experimental data (Materials and methods). Specifically, we rely on parameter estimation techniques for feedback systems to find a model that minimizes a quadratic error term between observed and predicted readouts, subject to simplicity constraints on the number of interactions in the system. Third, the fitted model can be used to predict the system's response to unseen perturbations (for example, combinations of drugs), and to gain new insight into the system's architecture. Testing modeling power for combinatorial perturbations in breast cancer cells Dual drug perturbation experiments in MCF7 breast cancer cells To directly test the power of the approach, we performed an independent experimental study in MCF7 human breast carcinoma cells. As perturbants of the system, we chose compounds targeting EGFR (ZD1839), mTOR (rapamycin), MEK (PD0325901), PKC-δ (rottlerin), PI3 kinase (LY294002) and IGF1R (A12 anti-IGF1R inhibitory antibody). As relevant readouts of molecular and phenotypic responses, we chose phospho-protein levels of seven regulators of survival, proliferation and protein synthesis (p-AKT-S473, p-ERK-T202/Y204, p-MEK-S217/S221, p-eIF4E-S209, p-c-RAF-S289/S296/S301, p-P70S6K-S371 and pS6-S235/S236) as well as flow cytometric observation of two phenotypic processes (cell cycle arrest and apoptosis) (Figure 2). Inhibitors were administered singly and in pairs, followed by EGF stimulation. When recording responses of protein phosphorylation, we used the average response at 5 and 30 min as the surrogate for steady-state values. To build models, we represented the state of each of the above perturbation targets (signaling proteins), as well as each of the readouts, by one state variable yi. We then used the proposed optimization procedure (Materials and methods) to estimate the coupling parameters wij and other parameters, resulting in predictive models of response in terms of these system variables. Figure 2.Breast cancer cells as a multiple input–multiple output system. To generate data for model construction, we treated human MCF7 breast tumor cell lines with one natural ligand (epidermal growth factor (EGF)) and six inhibitors, singly and in combination. The treatment protocol used EGF, an IGF1 receptor inhibitory antibody (A12) and inhibitors of the signaling molecules EGFR, PI3K, mTOR, PKC-δ and MEK. The inhibitors are ZD1839, LY294002, rapamycin, rottlerin and PD0325901. In the perturbation matrix (top panel, columns=experiments, rows=perturbations), a blue box indicates the presence of a particular perturbation and white indicates absence. For each treatment, we used western blots to detect levels of the proteins phospho-AKT, phospho-ERK, phospho-MEK, phospho-eIF4E, phospho-c-RAF, phospho-p70S6K and phospho pS6. We used a FACS-based assay to quantify apoptosis (measured as the 'sub-G1 fraction,' Materials and methods) and G1 arrest (measured as the G1 fraction). Here, representative flow histograms depicting cell cycle distribution in MCF7 cultures treated with or without drug are shown (one experiment is shown, see Supplementary information for all measurements). Evaluation of predictive power: After model construction based on these experiments, we see good agreement between experimental observation of the response of the MCF7 cell line to the 21 different perturbations (top, columns=experiments, rows=readouts) and the model prediction (bottom). Statistical assessment is in Supplementary Table 1. For each readout, we quantify the system's response by the phenotypic index defined as log relative response in treated versus untreated cells. For some drug combinations, the phenotypic readout increases as a result of perturbation (orange), for others it decreases (blue). Download figure Download PowerPoint Quantitative prediction of system response We first assessed the predictive power of the derived models using leave-one-out cross-validation, in which one pair perturbation is left out of the analysis and then its effect predicted from information gained from all other perturbations. The resulting predictions were reasonably accurate for the nine different readouts. The best prediction was obtained for p-S6 phospho-protein levels (cross-validation error CV=0.02, Pearson correlation r=0.96) and the weakest for the G1 arrest phenotype (CV=0.07, r=0.45) (Figure 2 and Supplementary Table 1). We directly compared the performance of our modeling approach to one using a corresponding set of linear differential equations with the same optimization procedure. By comparison, predictions using the nonlinear approach agreed better with experimental observations for eight of the nine readouts. Using the nonlinear modeling approach, the prediction error was lower by up to 50% with correspondingly better correlation between predictions and experimental observations (Supplementary Table 1). Thus, we conclude that our method is capable of deriving reasonably accurate network models for the input–output behavior of MCF7 cells with respect to the readouts used. Detection of key regulatory mechanisms without prior knowledge From a set of perturbation experiments, how can one deduce the logical network structure of activating and inhibiting interactions between the key molecular components, similar to the familiar pathway diagrams in publications summarizing a set of molecular biological experiments? Here, we use the derived network models with the smallest global error (Etotal=ESSQ+λESTRUCT, Materials and methods) to infer causal connectivity diagrams. The inference is based on the assumption that interactions in sufficiently simple models that fit experimental observations, called 'good' models, represent an underlying causal relationship between system components modeled by the system variables yi. Such a relationship can be either an indirect regulatory effect or a direct physical interaction that would be observable in vitro with purified components. Using our Monte Carlo algorithm, we generated a population of 450 good models from the MCF7 dual drug perturbation experiments. From these, we assessed the statistical significance of the individual interactions both in terms of a posterior probability (which is obtained directly from the Monte Carlo process, see Materials and methods) and a 90% confidence interval constructed by boot-strapping simulations (Table I). We now discuss the connectivity of the best model, i.e. the one with the smallest error (schema in Figure 3, explicit equations in Materials and methods) relative to the known biology of regulatory pathways in the MCF7 breast cancer cell line. Figure 3.Use of MIMO models to infer regulatory interactions in breast cancer cells. The interaction matrix wij from a set of good models can be used to infer regulatory interactions (squares=inputs; circles=internal system variables and other observables). Positive wij means activation and negative wij means inhibition of the target. Interestingly, some of the interactions are well known in MCF7 cells (green arcs) and others constitute predictions (orange arcs). See Table I for functional comments on interactions. No underlying pathway model was used—the network is a straightforward interpretation of the optimized model parameters wij. The EGFR → MEK → ERK and PI3K → AKT, canonical pathways are identified. Also, note the detection of self-inhibitory interactions in MEK/ERK signaling, identification of eIF4E and AKT as direct regulators of apoptosis and G1 arrest. Download figure Download PowerPoint Table 1. Statistical assessment of inferred interactions in MCF7 cells Regulator Target Mean interaction 90% CI Posterior probability (%) Comment EGFR → MEK 0.78 (0.67, 0.88) 99 0 Activation of MEK by RTK signaling MEK → p-ERK 0.39 (0.26, 0.48) 99 0 MAPK signaling PI3K → p-AKT 0.57 (0.47, 0.69) 93 0 PI3K-dependent activation of AKT PKC-δ → p-RAF 0.25 (0.17, 0.32) 71 0 Prediction mTOR → p-P70S6K 0.25 (0.20, 0.29) 41 0 mTOR signaling pathway p-ERK → p-eIF4E 0.19 (0.14, 0.25) 38 0 Prediction p-ERK → p-ERK 0.29 (0.22, 0.40) 35 0 Prediction p-P70S6K → p-S6 0.54 (0.53, 0.56) 20 0 S6 kinase activates S6 Apoptosis → G1 arrest 0.37 (0.30, 0.44) 20 0 Prediction G1 arrest → EGFR 0.39 (0.29, 0.49) 13 0 Low significance prediction MEK → p-RAF 0.29 (0.16, 0.41) 13 0 Low significance prediction p-eIF4E → PI3K 0.05 (0.02, 0.08) 13 1 Low significance prediction p-AKT Apoptosis −0.39 (−0.43, −0.36) 0 77 AKT controls survival in MCF7 cells PKC-δ G1 arrest −0.12 (−0.20, −0.06) 0 41 Prediction p-ERK EGFR −0.28 (−0.42, −0.17) 0 36 MAPK negative feedback loop G1 arrest p-eIF4E −0.33 (−0.36, −0.30) 0 32 Predicted bi-stable regulation p-eIF4E G1 arrest −0.19 (−0.25, −0.12) 0 20 Predicted bi-stable regulation Apoptosis p-P70S6K −0.39 (−0.42, −0.36) 0 19 Low significance prediction IGF1R G1 arrest −0.06 (−0.14, −0.01) 0 16 Low significance prediction G1 arrest p-ERK −0.09 (−0.16, −0.01) 1 15 Low significance prediction mTOR Apoptosis −0.07 (−0.12, −0.02) 1 12 Low significance prediction MEK MEK −0.13 (−0.21, −0.07) 0 8 Low significance prediction EGFR p-RAF −0.03 (−0.11, −0.04) 6 2 Low significance prediction Statistical assessment of inferred molecular interactions as shown in Figure 3. Column 'Mean Wij' shows the average interaction strength between the target and the regulator, as obtained from 200 bootstrapping simulations (see Supplementary information). 90% confidence intervals (CI) for the interaction strength are shown. The left column of posterior probabilities shows the fraction of models with a positive regulation in the Monte Carlo simulation. The right column shows the fraction of models with negative regulation (for instance, inhibition of apoptosis by p-AKT was present in 77% of models). The names p-P70S6K and p-p70S6K are synonymous. Interpretation of derived network structure In comparing the inferred connectivity with mechanisms known to occur in MCF7 cells (Table I), two caveats are important. (1) The logical nodes in our models are defined precisely as the perturbed and observed molecular species, i.e. the targets of drug perturbation and the targets of specific observed antibody reactions, and may not be exactly identical to a single molecular species. For example, 'EGFR' refers to the direct target(s) of activation by EGF and of inhibition by the drug ZD1839, and these two are assumed to be identical. (2) The models make no reference to unperturbed or unobserved nodes, e.g. whereas p-AKT is in the network model, the unphosphorylated AKT is not. With these caveats in mind, one can use the models both for confirmation and prediction of interactions. Of the 23 interactions in the best model, 14 had a posterior probability in the range of 20–99% (Table I). Of these, several statistically robust interactions clearly confirm canonical pathway structures. (i) The MAPK cascade downstream of the EGF receptor is detected as a chain of interactions between EGFR, MEK and ERK (Figure 3 and Table I). (ii) The negative feedback regulation of MAPK signaling is captured as negative interaction from ERK to EGFR, and as a moderately significant self-inhibition of MEK (see Discussion). (iii) PI3K-dependent signaling and the tendency for MCF7 cells to be dependent on AKT activation for survival are detected as interactions between PI3K, AKT and the apoptosis phenotype. (iv) The model inference that apoptosis is controlled by p-AKT, but not p-ERK, is in agreement with previous results in MCF7 cells (Simstein et al, 2003; DeFeo-Jones et al, 2005). (v) mTOR downstream signaling is detected as interactions between mTOR, p70S6K and ribosomal S6 protein (Mingo-Sion et al, 2005). The derivation of these expected interactions from a small set of perturbation experiments, without prior pathway knowledge, underscores the non-trivial value of the model building approach and provides some confidence in the concrete predictions of logical regulatory interactions for MCF7 cells (Table I), which are discussed below. Discussion In summary, our evaluation in breast cancer cells supports two main conclusions. First, our approach to model construction can be used to build reasonably accurate quantitative predictors of pathway responses to combinatorial drug perturbatio