Models from Experiments: Combinatorial Drug Perturbations of Cancer Cells
Sven Nelander,Weiqing Wang,Bengt‐Olof Nilsson,Qing‐Bai She,Christine A. Pratilas,Neal Rosen,Peter Gennemark,Chris Sander
DOI: https://doi.org/10.1038/msb.2008.53
IF: 13.068
2008-01-01
Molecular Systems Biology
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>at0 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