Network quantification of EGFR signaling unveils potential for targeted combination therapy

Abstract

The epidermal growth factor receptor (EGFR) signaling network is activated in most solid tumors, and small-molecule drugs targeting this network are increasingly available. However, often only specific combinations of inhibitors are effective. Therefore, the prediction of potent combinatorial treatments is a major challenge in targeted cancer therapy. In this study, we demonstrate how a model-based evaluation of signaling data can assist in finding the most suitable treatment combination. We generated a perturbation data set by monitoring the response of RAS/PI3K signaling to combined stimulations and inhibitions in a panel of colorectal cancer cell lines, which we analyzed using mathematical models. We detected that a negative feedback involving EGFR mediates strong cross talk from ERK to AKT. Consequently, when inhibiting MAPK, AKT activity is increased in an EGFR-dependent manner. Using the model, we predict that in contrast to single inhibition, combined inactivation of MEK and EGFR could inactivate both endpoints of RAS, ERK and AKT. We further could demonstrate that this combination blocked cell growth in BRAF- as well as KRAS-mutated tumor cells, which we confirmed using a xenograft model.

Figure 1

General outline of the study. A panel of six colon cancer cell lines was chosen, and profiled by sequencing selected cancer-related genes. The cells were then systematically perturbed with four kinase inhibitors and two ligands of growth receptors, and phosphorylation of key signaling proteins was measured using the Luminex platform. These data were used for parameterizing a mathematical model, which was then used to predict combinatorial treatments.

Figure 2

Generation of systematic perturbation data. (A) Perturbations consisted of two ligands (red nodes), and four pharmacological inhibitors (yellow flashes). These were applied alone and in inhibitor-ligand combinations for the indicated time points. Then eight phosphorylation signals were measured (blue nodes). (B) Log2 fold change of phosphorylation in response to the perturbations for the six indicated cell lines. Displayed response range was limited to ±3.5 (approx. 10-fold). Source data for this figure is available on the online supplementary information page.

Figure 3

EGFR modeling required network alteration that led to the discovery of a new effect of BMS345541 on ERK. (A) Scheme of the modeling pipeline from the starting network to its parameterization. Systematic perturbation data and a starting network serve as input. The core-fitting routine consists of three steps, and is illustrated by a four-node example network. First, non-identifiable parameters are detected and reparameterized. Second, the identifiable parameter combinations are fitted to the experimental data. Third, connections not significantly contributing to the fit are removed and subsequently the remaining parameters are refitted. The resulting network contains information about the strength, direction and sign of its connections. (B) Chi-squared values of the models trained on data of five cell lines. Dark bars indicate values for the initial model. Blue, red and yellow bars show values after model reduction of literature network, with a link from IGFR to RAF, or when additionally including a link from IKK to ERK, respectively. (C, D) Log2 fold change in time-series experiments after treatment with IKK inhibitors or solvent controls to investigate the IKK–ERK relationship in HCT116 cells. (C) IKK inhibitor BMS345541 treatment results in an increased phospho-ERK level, whereas treatment with IKK inhibitors PHA408, PS1145 and solvent control (DMSO) results in no increase. The inhibition of IkBa phosphorylation is comparable for all three inhibitors after 1 h (Luminex, n=1). (D) Phosphorylation of p70S6K, a cytoplasmic target of ERK, increases after 60 min of BMS345541 treatment (western blot, n=2) when compared with DMSO and PBS control. (E) Expression of ERK target genes EGR1 and FOS decreases after BMS345541 treatment (qrt–PCR, n=2). Source data for this figure is available on the online supplementary information page.

Figure 4

Model fit uncoveres differences in network topology and signaling strength. (A) Heat maps for the five cell lines showing log2 fold changes of phosphorylation (filled triangles) and the corresponding model simulation (open triangles). (B) Network structure derived from the model fit. Dashed lines indicate edges that are removed in those cell lines marked by filled circles (same colors and layout as in A). (C) Model-deduced parameter values representing signaling strength for identifiable paths, inhibitors and feedback/cross talk expressed as log2 response coefficients. Orange and blue boxes indicate feedback from ERK via EGFR to MEK, and cross talk from ERK via EGFR to AKT.

Figure 5

ERK–AKT cross talk is mediated by EGFR and independent of RAS or RAF mutation. (A) Increase of AKT phosphorylation in HCT116 (KRAS mutation) and HT29 (BRAF mutation) cells incubated with MEK inhibitor AZD6244 (0.1 μM) or its solvent control (DMSO) for 1 h before application of TGFa or BSA for 30 min. Effect can be also produced by 10 min EGF treatment in the presence of AZD6244 (1 μM). (B) Schematic view of the proposed ERK–AKT cross talk via the EGFR feedback (red), predicting no effect of MEK inhibition (flash) on stimulation of other growth receptors. (C) Increase of Akt phosphorylation compared with solvent control (DMSO) with combinatorial treatment of AZD6244 (5 μM) and different growth factors (10 min). (D) Response of phospho-ERK and AKT to 10 min EGF treatment for different preincubation times of AZD6244 (1 μM) in HCT116 cells. BSA and DMSO are solvent controls for EGF and AZD6244, respectively. All data measured using Luminex assays and shown as fold change; significant deviations are indicated. Brackets indicate significant one-sided t-test, P<0.05. Error bars indicate s.d. of n⩾3 samples. Source data for this figure is available on the online supplementary information page.

Figure 6

Model-derived combinatorial treatment options unveil effective combinations. (A) Model prediction of selected double inhibitions (red bars) and Luminex measurements (dark bars) for phosphorylation of AKT in HCT116 and HT29 cells in the presence of TGFa in log2 scale. Fits for single perturbations are shown on the left as blue bars. Error bars for modeled results represent s.d. of 100 noised simulations (see Supplementary Methods) (B) Cell proliferation (Xcelligence) of HCT116 and HT29 cells in response to single and double inhibition and control. Arrows mark time of treatment, and areas represent the maximal range of the growth curves, with dark lines reflecting the mean (n>=3). For panels A and B, concentration of AZD6244 was 1 μM. (C) Tumor growth of coloretal cancer cell DLD-1 cells transplanted into nude mice, which received a daily oral gavage of one of two different concentrations of MEK inhibitor GDC-0973, EGFR inhibitor erlotinib alone or combined. Error bars represent s.e.m. with n=10. Source data for this figure is available on the online supplementary information page.