Model-Based Design of a Decision Tree for Treating HER2+ Cancers Based on Genetic and Protein Biomarkers.

Human cancers are incredibly diverse with regard to molecular aberrations, dependence on oncogenic signaling pathways, and responses to pharmacological intervention. We wished to assess how cellular dependence on the canonical PI3K vs. MAPK pathways within HER2+ cancers affects responses to combinations of targeted therapies, and biomarkers predictive of their activity. Through an integrative analysis of mechanistic model simulations and in vitro cell line profiling, we designed a six-arm decision tree to stratify treatment of HER2+ cancers using combinations of targeted agents. Activating mutations in the PI3K and MAPK pathways (PIK3CA and KRAS), and expression of the HER3 ligand heregulin determined sensitivity to combinations of inhibitors against HER2 (lapatinib), HER3 (MM-111), AKT (MK-2206), and MEK (GSK-1120212; trametinib), in addition to the standard of care trastuzumab (Herceptin). The strategy used to identify effective combinations and predictive biomarkers in HER2-expressing tumors may be more broadly extendable to other human cancers.

Oncologists have long recognized that patients bearing tumors within similar histolopathologies and clinical features can respond vastly differently to chemotherapies. Data emanating from large-scale cancer genomics projects such as the Cancer Genome Atlas (TCGA) are now unveiling the molecular basis of this phenomenon. The genome of every cancer is essentially unique, such that even histologically identical tumors can emerge from completely nonoverlapping patterns of genetic alterations.1 Molecular information is thus increasingly used to inform treatment decisions, as oncology drug development has shifted from broadly cytotoxic therapies targeting DNA replication and metabolism, to molecularly targeted agents which inhibit specific molecular aberrations (oncogenes) and biochemical pathways. Meaningful responses to single targeted agents are, however, rare and often transient in nature.2 At the cellular level, robustness to drug intervention arises from the multiple redundancies and feedback regulatory circuits embedded in the oncogenic pathways driving cancer cell growth and survival.3 Realizing the full potential of targeted agents thus requires approaches to rationally combine therapies and select patients most likely to respond to such combinations.4 Quantitative models relating the molecular features of tumors to pharmacologic response patterns are thus essential for achieving the vision of precision medicine.

Mechanism-based models of signaling transduction, often based on mass action kinetics-based ordinary differential equations (ODEs), have proven extremely valuable in advancing our fundamental understanding of cell biology.5 Such models have also enabled the rational design of drugs that inhibit aberrant receptor signaling in cancer.6–9 However, the output is typically short-term biochemical changes, limiting the ability to simulate more clinically relevant tumor responses which occur over days to weeks. On the other hand, empirical models of tumor growth kinetics have been employed for decades as tools to interpret dose–response relationships and compare agents, and more recently for predicting outcomes of clinical trials based on preliminary tumor size measurements.10,11 Such empirical formulations, however, are incapable of predicting drug combination effects or response biomarkers, two pressing needs in oncology. To address such limitations, hybrid translational models have more recently been developed to bridge the gaps between these two classes, using established knowledge of molecular cell biology to link pharmacodynamic responses to cell proliferation and survival.12,13 These have proven capable of extrapolating clinically effective dosing regimens and pharmacodynamic biomarkers from preclinical data,14–19 and predicting combination effects between novel agents.20,21 A common limitation of such models arises from the use of a single, or at most a small number of immortalized cancer cell lines to serve as prototypic experimental models of the disease, neglecting the extensive patient-to-patient variability in the molecular features of tumors and responses to therapy.

The ErbB network is one of the most extensively studied areas of signal transduction, serving as an archetype for both systems modeling and molecular oncology.22 The ErbB family proteins consist of four receptors (ErbB1–4, also known as HER1–4). In response to extracellular ligand binding, combinatorial patterns of receptor dimerization result in activation of the phosphatidylinositol-4,5-bisphosphate 3-kinase/Akt (PI3K) and mitogen-activated protein kinase/Erk (MAPK) signaling cascades, which promote cell survival and proliferation.23 ErbB signaling is frequently dysregulated in human cancers, one of the most common mechanisms of which is amplification of the ERBB2 gene,24 resulting in overexpression of the surface receptor HER2 and constitutive signaling. The monoclonal antibody trastuzumab (Herceptin; Genentech/Roche, South San Francisco, CA) is standard of care in HER2+ disease, and could be considered the prototype for molecular targeted drugs. Many HER2+ cancer patients, however, do not respond to Herceptin, and those that do respond initially often develop resistance and relapse, particularly in the metastatic setting.25 Exposure to the HER3 ligand heregulin (HRG; NRG1), for example, induces HER2/HER3 heterodimer signaling and resistance to Herceptin in preclinical models.26 A wide variety of targeted therapies are currently undergoing clinical evaluation for treating Herceptin-refractory HER2+ disease, including tyrosine kinase inhibitors of HER2, PI3K, and MAPK pathway components, and biologics against HER3.27 However, the molecular and genetic determinants of sensitivity to such combinations remain obscure.

We have previously developed a semi-mechanistic model of ErbB signaling in HER2+ cancer cells, linking receptor engagement, via PI3K and MAPK pathway activation, to tumor growth kinetics. This was, however, parameterized using data from a single breast cancer cell line (BT-474), which is largely dependent on PI3K signaling. While the PI3K and MAPK pathways are aberrantly activated in the majority of human cancers,28 addiction to either (or both) varies widely both across and within cancer indications.29 We wished to explore whether differential dependence on PI3K vs. MAPK signaling cascades affect pharmacologic response patterns, and whether this information could be used to inform treatment strategies for HER2+ disease.

Based on these objectives we devised the following simulation-based strategy (summarized in Figure 1). Starting with our previously published model (representing a single prototypic cell line), we created two functional classes of cancer cells, PI3K or MAPK pathway-dependent, by tuning a logic gate linking the activation of these cascades to cell growth. Within each group, model parameters representing protein expression, gene mutations, and biochemical rates were randomized so as to generate synthetically heterogeneous populations. Randomized parameters thus serve as putative molecular biomarkers. Responses to 32 combinations of five clinically relevant drugs were then simulated across the synthetic populations, and optimal treatments defined based on median tumor growth inhibition. The PI3K and MAPK groups were then subdivided based on response (or lack thereof) to the population-optimal treatments. Biomarkers defining the outlier nonresponding cells were identified by comparing the underlying model parameters to that of the responsive group, and the biomarker-defined resistant groups subsequently analyzed to identify active drug treatments. Model predictions were assessed using in vitro proliferation assays in panels molecularly characterized cell lines, and clinical relevance using publicly available cancer genomic data.

Figure 1

Simulation-based strategy to identify effective drug combinations and predictive biomarkers using a mechanistic model.

Through this strategy, we defined molecular subtypes within HER2+ cancers based on a combination of genetic mutations (PIK3CA and KRAS) and microenvironment context (the HER3-ligand heregulin). By matching optimal drug combinations to these molecular subtypes, we specify a six-arm decision tree to stratify patients for treatment with combinations of the agents MM-111, lapatinib, MK2206, and GSK1120212. While this study was focused on five specific drugs within HER2+ cancer, we believe the approach may be more widely applicable to designing treatment strategies in other cancer indications.

RESULTS

Triple targeting of the HER2/HER3 complex is broadly effective regardless of downstream pathway dependence

Human cancer cells display a diverse spectrum of dependencies on PI3K and MAPK signaling pathways, and thereby sensitivity to agents targeting canonical components of these cascades, their upstream regulators, and downstream effectors. Here we sought to first explore whether dependencies on these two pathways in HER2-overexpressing cancers predict sensitivity towards combinations of targeted inhibitors. Five such agents were considered; the standard of care Herceptin, the HER2-targeted small molecule tyrosine kinase inhibitor (TKI) lapatinib (Tykerb, GSK), an antibody-based HER3 inhibitor MM-111,7 and small molecule inhibitors against the canonical PI3K and MAPK cascades MK2206 (an AKT inhibitor; AKTi), and GSK1120212 (trametinib, an MEK inhibitor; MEKi). A previously published computational model connecting HER2–HER3 signaling, via PI3K/AKT and MAPK/ERK signaling cascades, to tumor growth was utilized to assess all 32 possible combinations of the five agents.21 Two distinct classes of cancers were simulated based on preferential dependence of cell growth on the PI3K vs. MAPK cascades. These functional classes were specified by tuning parameters in the logic gate connecting pathway activity to cell growth regulation. Within each of the two classes, interindividual variability was simulated by randomizing 14 parameters (protein and gene-based putative biomarkers) within biologically feasible ranges (). Monte Carlo simulations were then used to assess growth inhibitory responses to the drugs across the population, and combinations rank ordered by median antitumor efficacy (Figure 2).

Figure 2

Triple targeting of the HER2/HER3 complex with MM-111, lapatinib, and Herceptin is broadly effective in both PI3K and MAPK-dependent cells. In silico screening of relative tumor growth in PI3K vs. MAPK-dependent cells in response to 32 combinations of five drugs (black = present, white = absent from regimens). MM-111, lapatinib, Herceptin (MHL) combination is indicated by asterisk. Relative growth (1 – Tumor growth inhibition).

Unsurprisingly, the dominant pattern that emerges is that layering on more drugs is predicted to increase efficacy. The model does not account for toxicity, which places constraints on the number of therapies, the specific combinations, and doses administered. While we cannot simulate toxicity associated with the various regimens, a few rules of thumb can be gleaned from clinical experience to date with combinations of targeted inhibitors. In general, combinations of biologics rarely produce unexpected or off-target toxicities, while combinations of small molecule kinase inhibitors often produce enhanced single agent-associated and unexpected toxicities.30

Examining the contribution of individual agents, as expected the AKTi-containing regimens were significantly more effective in PI3K-dependent tumors, while MEKi-containing regimens were more effective in the MAPK-dependent tumors. Regardless of pathway dependence, triple targeting the HER2/HER3 complex with the combination of MM-111, Herceptin, and lapatinib (MHL) was the most effective three-drug regimen when considering median growth inhibitory responses. We have previously assessed the tolerability of the MHL combination in mice as compared to other active regimens (AKT + MEK inhibition), and found it to be both well tolerated and effective preclinically (using a BT-474 xenograft model of HER2+ breast cancer), producing synergistic antitumor activity as compared to single agents.21 Twenty-nine patients with HER2+ solid tumors have also been treated with the regimen as part of a multiarm phase I safety study. Adverse events were similar to standard of care, and no maximum tolerated MM-111 dose was identified.31

While the MHL combination appears to be well tolerated and effective in both classes of HER2+ cancers, large variability in response to the regimen (3 orders of magnitude) was observed within both synthetic populations. To identify biomarkers underlying this response variation, in each class we separated nonresponding outliers (top 10%) from the median responding population (bottom 50%), and examined differences in parameter values (biomarkers) between the groups. Heightened expression of the ligand heregulin (HRG) emerged as a resistance mechanism shared by both tumor subclasses, consistent with its established role in mediating both Herceptin and lapatinib resistance26 and the mechanism of action of MM-111 as a competitive ligand antagonist. The top predictors of resistance in the PI3K and MAPK-dependent cells were found to be constitutive activating mutations within the respective signaling cascades (P = 2 × 10−16 and 3 × 10−15, rank-sum test; Figure 3). It is notable that biomarkers of resistance to MM-111, lapatinib, and Herceptin monotherapies were quite diverse, and the MHL combination markers were not linear derivatives from the individual agents ().

Figure 3

Biomarkers predictive of MHL regimen resistance in PI3K and MAPK-dependent tumors. (a) Cartoon schematic of ErbB signaling network model, connecting receptor activity, via PI3K and MAPK cascades to cell growth regulation. The 14 model parameters selected to vary in the Monte Carlo simulations are indicated numerically. Blue arrows indicate stimulatory relationships, and red inhibitory. (b) Median differences between each of the 14 model parameters (“biomarkers”) between the top 10% of resistant vs. 50% responsive tumors were compared. Shown are the Rank Sum P-values of those differences, expressed as –Log10. PI3K and MAPK pathway activating mutations emerge as the most significant biomarkers in the two classes.

PI3K and MAPK pathway activation mutations frequently co-occur in HER2+ disease

To assess the clinical relevance of these genetic alterations in HER2+ cancer, copy number and mutational profiles for the genes ERBB2, and commonly mutated components of the PI3K (PIK3CA, PIK3R1, and PTEN) and MAPK cascades (BRAF, KRAS, HRAS, NRAS) were extracted from primary tumor genome sequence data in the Cancer Genome Atlas (TCGA) using cBioPortal.32 We then examined co-occurrence of HER2-amplification with nonsynonymous PI3K and MAPK pathway mutations (Figure 4a). The majority of primary HER2+ cancers harbored significant frequencies of PI3K and/or MAPK pathway mutations, relatively consistent with frequencies of de novo Herceptin resistance.26 Samples in TCGA are largely primary tumor resections, and thus treatment-naive. As activating mutations in PI3K and MAPK pathways are established mechanisms of acquired resistance to ErbB targeted therapies,33,34 we expect these frequencies to be significantly larger in Herceptin-refractory tumors.

Figure 4

Co-occurrence of ERBB2 gene amplification with PI3K and MAPK pathway activating mutations, and predicted treatment strategies for the respective subpopulations. (a) ERBB2-amplified cancers represented in TCGA were segregated by PI3K (PIK3CA, PIK3R1, PTEN) and MAPK (KRAS, HRAS, NRAS, BRAF) pathway activating mutations. Those with neither are annotated as “wildtype” (WT). (b) Simulated relative growth (1-TGI) of PI3K-mutant, MAPK-mutant, and wildtype (WT) cells to select combinations of the five drugs considered. Shown are median responses, error bars indicate population quartiles.

We next parsed the response profiles within the PI3K-mutant and MAPK-mutant groups to find combinations with greater activity in these genetically defined subpopulations. Simulations predicted that switching out lapatinib in the MHL combination for an AKT inhibitor or MEK inhibitor in the respective PI3K and MAPK-mutant tumors would be significantly more effective than the MHL triplet as compared to the wildtype counterparts (Figure 4b). We next sought to test these biomarker-response predictions experimentally using a panel of genetically characterized HER2+ cancer cell lines.

PI3K and MAPK-activating mutations confer resistance to HER2/HER3 inhibitor combinations, but sensitivity to AKT and MEK inhibitor-containing regimens

In vitro proliferative dose responses to combinations of MM-111, Herceptin, lapatinib, MK2206, and GSK112012 were assessed in seven HER2+ cell lines, in the presence and absence of 5nM exogenous HRG ligand. This panel included cells harboring PIK3CA mutations (E545K and H1047R), KRAS mutations (G12D/C), as well as PI3K and MAPK “wildtype” (WT) cells. One WT cell (NCI-N87) was genetically engineered to express the PIK3CA mutants, providing a syngeneic platform to assess the effect of these single point mutations. The E454K and H1047R mutants had identical effects on drug response patterns, and were thus treated as replicates in the analyses. The panel also covers a diversity of indications, including breast (BT-474, MDA-MB-175-VII, MDA-MB-361), stomach (NCI-N87), colon (COLO-678), and esophageal (KYSE-410) cancers.

Due to different growth kinetics, microenvironment conditions, and time frames, a linear relationship between in vitro (cell) and in vivo (tumor) growth responses to drug treatment is unlikely. However, if drug effects are cancer cell-autonomous (i.e., not dependent on microenvironment or immunologic effects) the rank order of treatment responses should be consistent. Thus, we may use relative in vitro growth inhibitory responses to assess model predictions. Consistent with predictions, both PIK3CA and KRAS mutations significantly reduced the activity of the MHL combination and correspondingly enhanced sensitivity to AKT and MEK inhibitors (Figure 5a). Focusing on the engineered NCI-N87 cells, the single point mutations were sufficient to shift sensitivity from the lapatinib to the AKTi-containing treatments. Also note that sensitivity to AKT and MEK inhibition varied significantly between WT cells, representing differential dependence on the PI3K and MAPK pathways.

Figure 5

Resistance to PI3K and MAPK activating mutations can be overcome by switching lapatinib for AKT or MEK inhibitor-containing regimens, respectively. (a) In vitro proliferative responses (% growth inhibition at 1uM) to select drug combinations across seven cell lines, with and without stimulation by 5 nM heregulin (HRG). The heatmap is organized and labeled by treatment regimens (lapatinib (LAP), MK2206 (AKTi), and GSK1120212 (MEKi)-containing) and cell genotype (WT, KRAS, PIK3CA-mutant). (b) Average regression coefficients (BETA) for each drug across the three genetic ally defined cell populations, +/− heregulin (HRG). (c) Treatment decision tree based on the above results.

To distill the drug–response relationships from this multidimensional data (only a subset of treatments at a single dose (1 µM) is shown in Figure 5a; full data are provided in Supplemental File S1), multivariate regression models were parameterized for each cell line, simulating percent cell growth inhibition as linear combinations of the five input drug concentrations. Regression coefficients (BETAs) thus capture the relative responsiveness of each cell to the individual drugs. Average values of these five parameters were computed for each genetic subgroup (WT, KRAS-, and PIK3CA-mutant) in the presence and absence of exogenous HRG stimulation, resulting in 5 × 2 × 3 = 30 summary parameters (Figure 5b). This model formalism was chosen due to its simplicity, as a data compression tool (7 cell lines × 124 treatments × 4 replicates = 3,472 measurements, reduced to 30 summary parameters). The model implicitly assumes drug effects display log-linear dose responses and additive effects. Examination of the residuals reveals these assumptions are not drastically violated (). The model parameters (BETAs) can thus be used to evaluate consistency between model predictions and the experimental data.

Across all three genetic subgroups HRG stimulation increases sensitivity to MM-111, and decreases sensitivity to all four remaining drugs, particularly lapatinib. This is consistent with the known role of HRG-HER3 signaling as mediating resistance to multiple anticancer therapies (including HER2 inhibition),35 and the mechanism of action of MM-111 as a competitive HRG antagonist.7 Consistent with model predictions (Figure 4b), KRAS-mutant cells are the most sensitive to the MEK inhibitor, WT cells to lapatinib, and PIK3CA-mutant cells to the AKT inhibitor. The differential sensitivity of the PIK3CA mutant cells to AKT inhibition, however, is less pronounced than predicted. We attribute this to a limitation in our model's simplifying assumption as AKT being the sole effector of PI3K signaling. While AKT is the canonical output of the cascade, PI3K also activates MAPK signaling36 and MTOR.37 It is also notable that despite its well-established clinical benefit, Herceptin showed relatively little in vitro activity in any of the cell lines tested. Herceptin's activity in vivo is in part attributable to immunological rather than cell signaling effects,25 a facet not captured using in vitro proliferation assays.

We devised a simple decision tree summarizing our results as a guide for treating HER2+ disease with these agents (Figure 5c). As Herceptin (rather than lapatinib) is standard of care in HER2+ cancers and well tolerated, this is a backbone of all combination regimens. Relative sensitivity to the three small molecule kinase inhibitors (KI) considered was determined by genetic status of the PI3K and MAPK pathways, and due to the expected toxicity of small molecule KI combinations,30 we allowed for inclusion of only one of the three as part of a combination. Regardless of genetic status, MM-111 overcomes HRG-mediated therapy resistance, and should therefore be included if HRG is present. It is notable that this treatment stratification is based on a combination of both cell-autonomous genetic (PI3K and MAPK mutations) and microenvironment (HRG) determinants.

The accuracy of the decision tree in assigning maximally effective treatment (see File S1) for each of the 14 samples considered in Figure 5a is 50% (compared to 1/6 expected by chance, P = 4.1 × 10−3; binomial test), or 93% in choosing one of the top two most effective treatments (13/14 vs. 1/3 expected by chance, P = 6.1 × 10−6). Misclassifications are largely attributable to an added benefit of MM-111 in the absence of exogenously provided HRG, which may emanate from autocrine secretion of the growth factor, and/or inhibition of ligand-independent HER2/HER3 signaling. Taking this into account, accuracy in assigning optimal treatments is 79% (11/14 vs. 1/3 expected by chance, P = 6.9 × 10−4).

DISCUSSION

Through the strategy outlined in Figure 1 we were able to define pharmacologically distinct subtypes within HER2+ cancers based on a combination of genetic mutations (PIK3CA and KRAS) and microenvironmental context (the HER3-ligand Heregulin). Optimal drug combinations were matched to these molecular subtypes in silico, validated in vitro, and the results used to design a six-arm decision tree for the treatment of HER2+ cancers with combinations of MM-111, lapatinib, MK2206, and GSK1120212. While the focus of this study was limited to HER2+ cancers and the five drugs considered, the general framework and approach may be expandable to other cancer indications. Success largely depends on choosing appropriate and representative experimental platforms; the cell lines, pathway inhibitors, and quantitative molecular assays necessary for model parameterization and testing. Given the range of modeling approaches available to draw from38 the “optimal” choice, particularly related to the degree of mechanistic detail to include and modeling formalism, depends on the specific questions at hand and data available.39,40 Given the large gaps in current understanding of cell signaling networks,41 we bias toward using minimal mechanistic details and embedded assumptions, and gradually built out complexity as dictated by necessity.

Key predictions emanating from the model simulations were validated by in vitro pharmacologic response patterns. First, the ErbB3 ligand heregulin (HRG) drives resistance to HER2-inhibition, which can be overcome with an HRG antagonist such as MM-111. Second, constitutive activating mutations in PIK3CA and KRAS can drive resistance to dual HER2/HER3-targeted therapies, which can be overcome by switching the HER2 inhibitor lapatinib for AKT or MEK inhibitors. Genomic profiles from primary tumor samples suggest these genetic mechanisms to be clinically relevant within HER2+ disease. While not tested experimentally in this study, tumors harboring activating mutations in both cascades may require treatment with both AKT and MEK inhibitors, a possibly effective yet poorly tolerated combination.42

The use of individually tailored drug combinations has been widely touted as a solution to improve the effectiveness of anticancer drug therapy.43 However, the design of such biomarker-based regimens through empirical screening alone is largely unfeasible. Consider that the Cancer Cell Line Encyclopedia (CCLE) annotates 57 cell lines as ERBB2-amplified (copy number ≥ 4N).44 Molecular information available on these cells includes targeted sequencing of 1,060 cancer-associated genes, and whole-genome mRNA expression profiles (∼20,000 genes). This exemplifies the “large P small N problem”; the number of features (potential biomarkers) far exceeds the number of samples available for testing. Paradoxically, the availability of Big Data makes data-driven computational approaches alone insufficient for biomarker identification. Quantitative systems pharmacology models provide an alternative platform for rapidly screening hypotheses in silico, and thus focusing efforts on the most relevant experiments to perform. Model simulations guided our selection of the specific drug combinations, and the genetically defined cell lines to test the effect of PIK3CA and KRAS mutations on drug combination response patterns.

The decision tree (Figure 5c) neatly summarizes our experimental results, but ultimately may be useful to guide clinical development of the combination regimens. Patients are diagnosed with HER2+ cancer based on the immunohistochemistry-based HercepTest and FISH-based gene amplification assays performed on biopsied tumor tissue. Multiple additional assays may be run on the tissue sample, including an assessment for mutations in the genes PIK3CA and KRAS (or other pathway activating mutations such as in PTEN or BRAF), and expression of HRG (either protein via immunohistochemistry, or the surrogate mRNA transcript NRG1 via quantitative reverse-transcription polymerase chain reaction (qRT-PCR) or in situ-hybridization). Based on the results of such tests, patients would be classified into the six HER2+ subcategories defined in Figure 5c. As genomic profiling is becoming an increasingly routine part of clinical practice45 and NRG1 transcript expression can be quantified from tumor biopsies,46 measurement of the three biomarkers is indeed clinically tractable. While measuring the biomarkers in tumor tissue may be feasible, clinical trial designs capable of assessing multiple hypotheses (the efficacy of six combination treatments, and the predictive value of three biomarkers) would be required to clinically validate our findings. The optimal design of such multiagent, multibiomarker “basket trials” remains an open problem in clinical oncology,47 well beyond the scope of this study. So while multiple factors are required to advance preclinical results such as these into clinical practice, we believe that our work serves as proof of principle that model simulations, incorporating clinically relevant genetic and protein biomarkers, can be employed to advance the vision of precision medicine in oncology.

METHODS

Translational model and in silico drug combination screening

The semimechanistic model connecting ErbB receptor signaling, through PI3K/AKT and MAPK/ERK cascades to tumor growth is depicted schematically in Figure 3a. The underlying assumptions, mathematics, and data used for model parameterization have been previously described in extensive detail21 and are summarized in the Supplementary Methods. A logic OR gate was implemented to specify cell survival as dependent on phopsho-AKT and -ERK, the relative balance determined by empirical weighting parameters w and w. PI3K pathway dependence was thus simulated by setting w = 0.95 and w = 0.05, and conversely for MAPK pathway dependence. Tumor heterogeneity was simulated via Monte Carlo sampling of the model parameters (prospective biomarkers) listed in from log-uniform distributions.

For drug screening, tumor growth was simulated over 2-week periods, with lapatinib administered daily at 1250 mg and MM-111 weekly at 20 mg/kg, as per clinical regimens. As MK2206,48 GSK1120212,49 and Herceptin50 are cleared relatively slowly, and given uncertainties in their PK-PD relationships, the pharmacokinetics of these agents were ignored and in vivo effects simulated by constant target suppression (95% inhibition of phospho-AKT, -ERK, and -HER2). Tumor growth inhibition (TGI) is defined as the change in tumor size from baseline 2 weeks posttreatment as compared to matched untreated control. The 2-week timepoint chosen for analysis is shorter than the typical time between tumor scans (bimonthly); however, longer simulated treatment durations did not affect rank order of treatment effects, the metric used to experimentally assess predictions.

In vitro cell growth assays

Cellular responses to MM-111, Herceptin, lapatinib, MK-2206, and GSK1120212 were evaluated by CellTiter-Glo luminescent cell viability assays (Promega, Madison, WI). Cells were seeded at 700 cells per 384-well plate in 10% fetal bovine serum (FBS) cell growth medium and treated with the five drugs separately and in combination at 1, 0.1, 0.01, and 0.01 μM, with and without 5 nM HRG-b1 prestimulation (for 4 hours). Cell viability was determined 72 hours posttreatment, and Cell Growth Inhibition (CGI) defined as cell viability in the treated conditions compared to matched, untreated control.

Multivariate linear regression models of drug combination effects

CGI was described using a multivariate linear regression function of the Log10 drug concentrations (C): where N = number of input drugs (5: MM-111, lapatinib, Herceptin, MK-2206, and GSK-1120212), and β = regression coefficients, estimated by maximum likelihood estimation, and errors by computing the Hessian matrix. Goodness-of-fit was assessed by examining the distribution of residuals and their correlations with input parameters ().

All model simulations and computational analyses were performed in MATLAB R2013b (MathWorks, Natick, MA).