Introduction

Synergistic drug combinations have been extensively explored for enhanced therapeutic efficacies1,2,3,4,5,6,7,8,9. In discovering and investigating synergistic drug combinations, the level of synergism is typically measured and quantified by the drug combination index (CI, a quantitative measure of drug combination effects defined in Method Section) such as Chou and Talalay’s CI from experimental dose-response data1,3,10. Based on our literature search study, over 523 papers since 2004 have reported the discovery and optimization of synergistic drug combinations based on the experimentally determined CIs. In-silico tools that can predict CIs without the time-consuming and costly measurement of dose-response data are highly useful for facilitating the discovery of synergistic drug combinations.

Computational methods have been developed for predicting drug combination effects from gene expression profiles of drug-treated samples11,12,13,14,15 and simulation of drug-targeted signaling16,17,18,19,20,21 and metabolic22,23,24,25 pathways. In particular, simulation of drug-targeted pathways is potentially useful for predicting CIs17,26, as demonstrated by the successful applications of the chemogenomic profile based models27,28 and the statistically-inferenced network models29,30 for the prediction of synergistic effects of drug combinations. But the ability of the pathway simulation methods in predicting CIs has not been adequately tested against the observed values of multiple drug combinations targeting multiple target combinations. More tests are needed for determining what the existing mathematical models are capable of and what need to be further improved. These also provide useful knowledge for developing drug or drug combination targeted mathematical models for a number of pathways targeted by drugs and drug combinations (e.g. EGFR-ERK31,32,33,34,35, apoptosis36,37, NFκB16,17, Wnt19 and disease-relevant metabolic22,23,24,25 pathways).

In this work, we developed and tested a mathematical model of drug and drug combination targeted EGFR-ERK pathway (Fig. 1) based on the ordinary differential equation model of Hornberg38. The method for developing this model is provided in the method section. This pathway was selected for two reasons. First, several kinases in this pathway have been targeted by individual inhibitor drugs and drug combinations with available experimental drug response and CI data39,40,41,42,43,44. Secondly, it is one of the pathways with well-established mathematical models31,32,33,34,35,38, ideal for developing and testing drug-targeted pathway mathematical models. The kinase inhibitor drugs included in our mathematical model are EGFR, BRaf and MEK inhibitors, which together with their combinations have been clinically used or tested for the treatment of melanoma, colon, gastric, pancreatic, non-small-cell-lung-cancer (NSCLC) and other cancers39,40,41,42,43,44.

Figure 1: Drug-targeted EGFR-ERK pathway schema in this study.
figure 1

The EGFR, Raf and MEK inhibitor is represented by the small green, blue and yellow colored node with a letter D respectively.

The inhibitory effect of each drug against its target was measured by the percentage reduction of the integrated non-drug-bound target level at different drug concentrations (target dose response curve), and the concentration that induces 50% reduction was taken as the half maximal inhibitory concentration (IC50 value). The integrated non-drug-bound target level refers to the integral of the free target level over the first 2 hours of signaling stimulation. The anti-proliferative effect of each drug or drug combination was measured by the percentage reduction of the integrated phosphorylated ERK (ppERK) level (described below) with respect to the concentration(s) of the drug or drug combination (anti-proliferative dose effect curve)45, and the concentration(s) that induce 90% reduction was taken as the half maximal inhibition of growth (GI50 value) of the drug or drug combination (details in the Method Section). The integrated ppERK level refers to the integral of the ppERK level over the first 2 hours of signaling stimulation (the ppERK level typically returns to the basal level <2 hours after HGF stimulation46). The anti-proliferative effect of a drug or drug combination was measured by its induced reduction of the integrated ppERK level for the following reason: The level of cell proliferation is linked to DNA synthesis, which is positively associated with the level of the integrated ERK2 activity47, while the latter is correlated to the integrated ppERK level. However, a nonlinear correlation between the integrated ERK activity and cell proliferation rates is observed, with significant reduction of proliferation typically occurs when the ppERK level falls below 10% of its peak value48,49. Therefore, the concentration(s) of a drug or drug combination that induces 90% reduction of the integrated ppERK level was used to measure the GI50 value.

Our mathematical model was validated against the previously observed and simulated effects of the regulation of EGF, PP2A, MKP3 on ERK activities and the regulation of EGF on RasGTP, and its predicted anti-proliferate effects of EGFR, BRaf and MEK inhibitors were compared to the experimental cancer cell-line growth inhibition GI50 values. Additionally, the simulated efficiency of EGFR-MEK combination inhibitors was consistent with experimental results. It was then used to predict the effects of the combinations of EGFR-BRaf and BRaf-MEK inhibitors, and the computed anti-proliferate dose-effect curves were used to derive the drug-combination isobolograms and CIs based on Chou and Talalay’s formula10. The predicted CIs were compared to the experimental values for evaluating the performance of the model.

Results and Discussions

Model validation against the observed signaling dynamics and the reported simulation results in the absence of a drug

Our mathematical model, described in the Method section, was first validated by comparing the simulated signaling dynamics in the absence of drugs with the literature-reported observations and simulation results. We found that the simulated time-dependent protein concentration profiles are in reasonable agreement with the experimentally determined profiles and the reported simulation results. At 8 nM EGF, the simulated ERK activation peaks within 5 minutes and decays within 2 hours (Supplementary Figure S1), which is consistent with the observed46 and the simulated38 ERK phosphorylation kinetics. The amount of simulated active RasGTP peaks at ~2 minutes and quickly decays within 10 minutes (Supplementary Figure S2), which is consistent with the observation that active RasGTP level in EGF- treated PC12 cells increases dramatically within 5 minutes and decays steeply within 10 minutes32. In simulating the regulation of EGFR-ERK signaling by phosphatases, we found that increasing PP2A level from 30000 to 60000 molecules/cell induces little changes in the peak amplitude but significantly reduces the duration of ERK activation (Supplementary Figure S3), while increasing MKP3 level from 9 × 106 to 1.5 × 107 molecules/cell significantly alters the peak amplitude and substantially reduces the duration of ERK activation (Supplementary Figure S4), which are consistent with the results of a reported simulation study showing that the duration of ERK activation is sensitive only to phosphatase reactions on MEK whereas the amplitude is most sensitive to phosphatase reactions on ERK50.

Model validation against the experimental anti-proliferation activities of individual drugs

Our mathematical model was further validated by the comparison of the predicted GI50 values of individual EGFR, BRaf and MEK inhibitor with the reportedly observed GI50 values of known inhibitors against different cancer cell-lines51,52. Figure 2 shows the computed anti-proliferative dose-response curves of the EGFR, BRaf and MEK inhibitor respectively, and the computed target inhibitory dose-effect curves of these three inhibitors are shown in Supplementary Figures S5–S7 respectively. From Fig. 2, the computed GI50 values for EGFR, BRaf and MEK inhibitor are 19, 113, 29 nM respectively, which are within 10 fold of the largest experimental GI values against the cancer cell-lines promoted by the EGFR pathway (45 and 53 nM for EGFR inhibitor gefitinib and erlotinib against non-small-cell lung cancer (NSCLC) cell-lines51, 1.1 μM for BRaf inhibitor sorafenib against NSCLC cell-lines50, and 14–200 nM for MEK inhibitor AZD6244 against thyroid cancer and melanoma cell lines52.

Figure 2: The computed anti-proliferative dose-response curve of individual inhibitors against EGFR-ERK pathway mediated cell growth signaling.
figure 2

(A) EGFR inhibitor; (B) BRaf inhibitor; (C) MEK inhibitor.

The Kd values of the EGFR, BRaf and MEK inhibitors were set such that it gives the median IC50 (10 nM, 30 nM, and 15 nM for EGFR, BRaf, MEK inhibitors respectively) corresponding to the highly potent IC50 values of the drugs frequently used in pharmaceutical research (12.5 ± 7.5 nM against EGFR, 28 ± 3 nM against BRaf, and 15.5 ± 3.5 nM against MEK, details in the Method Section). The predicted GI50 values of the EGFR and MEK inhibitors are substantially smaller than (by 4–6 fold) that of the BRaf inhibitor. This is consistent with the observed behavior of the highly potent inhibitors with the largest experimental GI values (by 5.5–22 fold) against the cancer cell-lines promoted by the EGFR pathway. Analysis of the interaction dynamics of our mathematical model (Fig. 1) suggested that this phenomenon arises from the additional bypass signaling through BRaf-CRaf dimers after BRaf inhibitor binding to BRaf. As a result, increasing activity of the BRaf inhibitor stimulates the bypass signaling mediated through CRaf, and vice versa, leading to an overall less efficient anti-proliferative effect by the BRaf inhibitor than that by the EGFR and MEK inhibitor respectively.

Our mathematical model, was derived from the conventional models31,32,33,34,35,38 and drug-target binding kinetics34 with the parameters fitted to the median IC50 values of potent inhibitors. Variation of the IC50 from 5 nM to 50 nM led to comparable degree of variations in the predicted GI50 values (9–95 nM, 19–195 nM and 10–99 nM for EGFR, BRaf and MEK inhibitor respectively). The consistency of the predicted and observed GI values suggests that our mathematical model may have some capability in facilitating the quantitative study of the anti-proliferative activities of the drugs and drug combinations targeting the EGFR pathway.

Model validation against the experimental anti-proliferation activities of EGFR-MEK inhibitors combination

We used our mathematical model to simulate the anti-proliferative effect and to compute the CI values of drug-pair EGFR-MEK inhibitors. As in the cases of drug combination studies1,3,10, the concentration-dependent anti-proliferative effect of drug-pair was measured by an isobologram that displays a curve of the dose-combinations of the two drugs that produce 90% reduction of the integrated ppERK levels within the first 2 hours of signaling stimulation. Figure 3 shows the computed isobologram for the combination of EGFR-MEK inhibitors. The computed CI values for the combinations of EGFR-MEK is 0.46. These are comparable to the experimental values 0.4–0.8 for EGFR-MEK40, which suggest that our mathematical model has some level of capability in facilitating the prediction of CI values.

Figure 3
figure 3

The computed isobologram for the combination of EGFR-MEK inhibitors.

Prediction of drug combination effects of EGFR-BRaf, BRaf-MEK

BRaf is a commonly mutated oncogene, yet there are no effective therapies exist for BRaf mutant cancers. Mostly because targeting BRaf or MEK kinase in BRaf mutant cancer, respectively, frequently activates other bypass pathways, thus constraining the effectiveness of these inhibitors as single-agents. In our mathematical model, we further predict the combination therapy strategies for BRaf mutant cancers with two drug combination effects, EGFR-BRaf and BRaf-MEK inhibitors, respectively. In modeling of these two drug-pairs, the same set of kinetic parameters for modeling as each individual drug was used. Figure 4 showed the computed isobolograms for the combination of EGFR-BRaf, BRaf-MEK inhibitors. The computed CI values for the combinations of EGFR-BRaf, and BRaf-MEK are 0.69 and 0.87, respectively.

Figure 4: The predicted isobologram for two drug combinations.
figure 4

(A) EGFR-BRaf inhibitors; (B) BRaf-MEK inhibitors.

The synergistic effects of these two drug combinations arise from the same type of synergistic mode of action: the complementary action involving positive regulation of a target by targeting two upstream-downstream points of a pathway8,53,54. For the EGFR-BRaf inhibitor combination, the inhibitory activity of the EGFR inhibitor reduces the activation BRaf to complement the inhibitory activity of the BRaf inhibitor. For the BRaf-MEK inhibitor combination, the inhibitory activity of BRaf reduces the activation of MEK to complement the inhibitory activity of MEK inhibitor. These two combinations have largely similar level of synergism as judged by the similar quantities of the computed and observed CI values, which is likely due to the same type of synergistic mode of action.

Prediction of combination effects of a farnesyltransferase inhibitor combined with a drug targeting BRaf or MEK

Mutated Ras genes that produce constitutively active Ras proteins are found to be involved in approximately 30% of human cancers, including several pancreatic and colon carcinomas. Anticancer therapeutic development has been focusing on blockage of hyper-activation of Ras protein, e.g. through farnesyltransferase inhibitors (FTIs) which inhibit the farnesylation of RAS, preventing Ras transform into the functional form that can be attached to the cell-membrane and subsequently transmit signals to the Raf-MEK-ERK (MAPK) signaling pathways, which have a major role in melanoma progression55. However, FTIs alone is frequently insufficient for the treatment of cancers. In our mathematical model, we simulated the combination effects of two types of drug combinations, FTI-MEK and FTI-BRaf inhibitors, respectively. FTIs are modeled in our study as drugs inhibiting the transformation of RasGDP to functional RasGTP. The association and dissociation constants kon and koff for FTIs are optimized to give IC50 of 10 nM, since most FTIs in pharmaceutical research has less than 10 nM IC50 values56,57.

Figure 5 showed the computed isobologram for the combination of FTI-MEK inhibitor, FTI-BRaf inhibitor. The computed CI values for the combinations of FTI-MEK inhibitor, and FTI-BRaf inhibitor are 0.74 and 0.78, respectively. The predicted synergism effects of these two drug combinations are consistent with recent literature report that combinations of lonafarnib (a specific FTI) with pan-RAF inhibitor sorafenib yielded additional growth-inhibiting effects than lonafarnib and sorafenib used alone58. Moreover, the apoptosis-inducing effect of BMS-214662 (another cytotoxic FTI) has been significantly enhanced through combination treatment with a MEK inhibitor, PD184352, leading to a complete suppression on invasive tumor growth in K562 and CD34 + CML cells59. Both our model and experimental studies suggest that the combination treatment of FTI and Raf and MEK-inhibitor may represent an effective alternative for melanoma treatment, and therefore worth further exploration.

Figure 5: The predicted isobologram for combination treatment of farnesyltransferase inhibitors (FTIs) and drugs.
figure 5

(A) FTI-MEK inhibitor; (B) FTI-BRaf inhibitor.

Conclusions

Our drug-targeted EGFR-ERK pathway mathematical model, developed based on existing drug-free models and drug-target interaction kinetics, and validated against a number of published experimental and simulation results, showed fairly good potential in predicting the individual and combination drug effects on the cell-proliferation signaling of the pathway. In particular, the predicted drug combination CI values are comparable to the experimental values for the combination of EGFR-BRaf, EGFR-MEK, BRaf-MEK, FTI-MEK, and FTI-Raf inhibitors41, suggesting that existing pathway models may be potentially extended for developing drug-targeted pathway mathematical models to predict drug combination CI values, isobolograms, and drug-response surfaces/curves as well as to analyze the dynamics of individual and combination of drugs. Our model provides a strategy to predict efficacy of these MAPK pathway inhibitors through novel targeted therapy combination approaches.

Further development efforts are needed to explore drug-targeted pathway mathematical models as practical tools for facilitating the discovery of synergistic drug combinations. For instance, some synergistic drug combinations are known to arise from the activities against multiple pathways8. Although several multi-pathway models have been developed60,61,62, there is a need to refine existing and develop more multi-pathway models to sufficiently cover the multiple pathways relevant to synergistic drug combinations. Moreover, in some cases, drug effects are influenced by pathway crosstalks63,64, and drug actions at the multiple sites65,66, states67, conformations1, and mutant forms1 of the same target. There is a need to add these elements into the existing and newly developed pathway models for extended coverage of therapeutically relevant biological networks and the modes of actions of individual drugs and drug combinations8.

Methods

Construction of Mathematical Model

Our pathway model schema is illustrated in Fig. 1. We based our model of EGFR signaling on that of Hornberg et al., which itself is a refinement of earlier work. The components of drug-target interactions were added to the ordinary differential equation model of Hornberg et al.38,68, with an additional ERK-CDC25C-EGFR feedback loop44 and the interactions between BRaf and CRaf69 incorporated. The model contains 126 distinct molecular species, and 190 elementary reactions; these reactions are described as a series of ordinary differential equations based on the mass action law. The model is parameterized by 123 kinetic parameters and 126 initial molecular concentrations. A typical enzyme catalyzed reaction in the pathway and the mass-action based rate is shown in the equations below, where the A and B represent species or concentration of species, depending on the context:

The rate equations are defined by the forward and reverse rate constants Kf and Kb or turnover Kcat used in the published models31,32,38,70,71,72,73 or reported from other literatures. The constituent molecular interactions, their kinetic constants, and molecular concentrations are detailed in Supplementary Table S1. Fourth order Runge-Kutta method with adaptive step-size control was used for solving these equations.

Kinetic Parameters

The types of parameters used in our mathematical model are protein-protein interactions, drug-protein interactions, and catalytic activities. The published simulation studies have shown that most parameters are robust and insensitive to significantly alter the overall pathway behavior31,32. Apart from the use of the parameters of the published mathematical models, additional parameters were obtained from the literatures based on the widely used assumption that the parameters measured in vitro and in some cell lines are generally applicable in most cases. For those protein-protein interactions with unavailable parameters, their parameters were putatively estimated from the known parameters of the relevant interacting domain profile pairs74,75 or other interacting protein pairs of similar sequences.

Drug targeting

Drugs targeting each component in the EGFR signaling pathway were simulated as forming complexes with target components based on the mass action law. These inhibitors include the individual EGFR, BRaf, MEK inhibitor or the combination of EGFR-MEK, EGFR-BRaf, and BRaf-MEK inhibitors respectively. For simulating the interaction between each drug and its target, we followed Bairy and Wong34 to introduce an extra species, the drug, and two reactions to describe drug-target association and dissociation, with the corresponding kinetic constants Kon and Koff determined such that the Koff takes a value of 0.01 s−1, while the Kon was determined such that the drug gives the median IC50 (10 nM, 30 nM, and 15 nM for EGFR, BRaf, MEK inhibitors respectively) corresponding to the highly potent IC50 values of the drugs frequently used in pharmaceutical research (Cetuximab 1 nM, Erlotinib 2 nM, Afatinib 14 nM and Gefitinib 33 nM against EGFR, Sorafenib 25 nM and Vemurafenib 31 nM against BRaf, and AZD6244 12 nM and RDEA119 19 nM against MEK).

Definition of IC50 and GI50 values

IC50 is defined as the half maximal inhibitory concentration that induces 50% reduction of integrated non-drug-bound target level. The integrated non-drug-bound target level refers to the integral of the free target level over the first 2 hours of signaling stimulation. For the reasons described in the Introduction Section, the drug concentration-dependent (from 0.001 nM to 10 μM) effects on the ERK-mediated cell-growth signaling process were determined by measuring the percentage of ppERK reduction within the first 2 hours of signaling stimulation. GI50 is defined as the growth inhibition concentration that induces 90% reduction of ppERK within the first 2 hours. Further analysis showed that the variation of the GI50 values from 80% to 99% of ppERK reduction resulted in insignificant changes in the predicted CI values (±0.1), and thus no changes in the qualitative conclusion of the synergistic or antagonistic levels of drug combinations. This indicates the robustness of our drug-binding models in predicting drug combination effects.

Computation of the combination index and isobolograms for quantitative determination of drug interactions

Quantifying drug interactions in drug combination studies and classifying the interactions into categories of synergy, additivity, or antagonism are of interest to many researchers. Isobologram and combination index (CI) analyses are widely used methods for evaluating drug interactions in combination cancer chemotherapy. The Loewe additivity model has been largely used as a reference model when the combined effect of two drugs is additive. The model can be written as in Equation (3):

where (D)1 and (D)2 are the respective combination doses of drug 1 and drug 2 that yield an effect of 50% growth inhibition, with (Dx)1 and (Dx)2 being the corresponding single doses for drug 1 and drug 2 that result in the same effect, which is by definition the GI50 of drug 1 and drug 2. When Eq. 3 holds, it can be concluded that the combined effect of the two drugs is additive. Based on Eq. 3, the combination index, defined in Eq. 4, can be used to classify drug interactions as synergistic, additive, or antagonistic.

A CI of less than, equal to, and more than 1 indicates synergy, additivity, and antagonism, respectively.

We simulated isobolograms for a pair of drugs with eight equally effective dose combinations for a particular effect level of GI50. GI50 normalized doses of drug 1 and drug 2 that give this effect in combination are plotted as axial points in the isobologram graphs. According to Eq. 4, the isobologram curves are expected to be parallel to the diagonal for additive drug pairs, concave for synergistic drug pairs, and convex for antagonistic drug pairs. We mainly concerned about the qualitative shape of the isobolograms for correctly identifying the drug pair category, and use the smallest CI of the eight drug dose combinations as the CI for this drug pair.

Additional Information

How to cite this article: Huang, L. et al. Predicting Drug Combination Index and Simulating the Network-Regulation Dynamics by Mathematical Modeling of Drug-Targeted EGFR-ERK Signaling Pathway. Sci. Rep. 7, 40752; doi: 10.1038/srep40752 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.