|Home | About | Journals | Submit | Contact Us | Français|
Targeted therapies have shown significant patient benefit in about 5–10% of solid tumors that are addicted to a single oncogene. Here, we explore the idea of ligand addiction as a driver of tumor growth. High ligand levels in tumors have been shown to be associated with impaired patient survival, but targeted therapies have not yet shown great benefit in unselected patient populations. Using an approach of applying Bagged Decision Trees (BDT) to high-dimensional signaling features derived from a computational model, we can predict ligand dependent proliferation across a set of 58 cell lines. This mechanistic, multi-pathway model that features receptor heterodimerization, was trained on seven cancer cell lines and can predict signaling across two independent cell lines by adjusting only the receptor expression levels for each cell line. Interestingly, for patient samples the predicted tumor growth response correlates with high growth factor expression in the tumor microenvironment, which argues for a co-evolution of both factors in vivo.
The combination of Herceptin® with chemotherapy demonstrated a dramatically increased survival benefit for a subset of women with HER2 amplified advanced breast cancer, which ultimately led to FDA approval in 1998.1 Since then, targeted cancer therapies have become an accepted therapeutic modality for the treatment of cancer and have contributed to a decrease in cancer related mortality.2 However, the benefit of targeted therapies to date has been restricted to 5–10% of solid tumors addicted to oncogenes.3–5 Identifying these relatively rare patients via predictive diagnostic tests relying on genomic biomarkers has created Precision Medicine.6–8
Retrospective analyses of several clinical studies of breast, gastric or lung adenocarcinoma identified increased receptor and/or growth factor expression as prognostic markers for patients with poor prognosis, which highlights the role of ligand-induced signaling as oncogenic drivers.9–12 Here we aim to decipher what drives ligand-induced proliferation.
We present the first comprehensive proliferation screen across 58 cell lines comparing to which extent the growth factors EGF (epidermal growth factor), HRG (heregulin), IGF-1 (insulin growth factor 1) and HGF (hepatocyte growth factor) induce cell proliferation. We find that about half of the cell lines do not respond to any of the ligands whereas the other half of the cell lines respond to a least one ligand. We compare the observed ligand-induced proliferation with the response to treatment with antibodies targeting the ErbB receptor family members, a subfamily of four closely related receptor tyrosine kinases (RTKs): EGFR (ErbB1), HER2/c-neu (ErbB2), HER3 (ErbB3) and HER4 (ErbB4) as well as the insulin growth factor receptor (IGF-1R) and the hepatocyte growth factor receptor (Met). Not surprisingly, the antibodies targeting the respective RTK inhibit ligand-induced proliferation. The antibodies also inhibited basal proliferation in some cell lines that do not respond to exogenous ligand addition, which could be driven by autocrine signaling.
The need has been recognized for computational approaches to deal with the complexity of signal transduction and its dysregulation in cancer to ultimately understand drug activity.13–17 Large collections of genetic and genomic data led to efforts to disentangle the complex mechanisms using machine-learning algorithms.18–21 It was previously shown that simulated patient-specific signaling responses derived from mechanistic signaling models using RNA sequencing data from patient biopsies can be robust biomarkers that are predictive of patient outcome.22 Here, we combined machine learning and mechanistic modeling to predict which cell lines proliferate in the presence of ligand. We used RNA sequencing data as inputs into a comprehensive mechanistic model capturing the ErbB, IGF-1R and Met signaling pathways. Our novel approach uses simulated signaling features and mutation status of a specific cell line as inputs into a Bagged Decision Tree, which predicts whether tumor cells proliferate in the presence of a growth factor. We achieved a substantial gain in accuracy compared to predictions based on RNA sequencing data alone by inclusion of simulated signaling features such as the area under curve of distinct heterodimers and phosphorylated S6 for in vitro models.
Applying this approach to patient data, the prediction of ligand-dependent tumor samples based on mRNA data from The Cancer Genome Atlas (TCGA) revealed that colorectal and lung cancer are the two indications most responsive to EGF, which agrees with the approval of EGFR inhibitors in these indications. In addition, the prediction of responders in patient samples revealed a correlation between predicted tumor growth and measured ligand expression in the tumor microenvironment, which argues for a co-evolution of ligand production and the ability of the tumor cells to respond to stimulation.
To investigate growth factor-induced proliferation we screened a panel of 58 cancer cell lines (10 ovarian cancer, 11 breast cancer, 13 lung cancer, 11 gastric cancer, and 23 colorectal cancer cell lines) for response to the exogenously added ligands EGF, HRG, HGF, and IGF-1 (Supplementary Fig. 1) that bind to EGFR, ErbB3, Met, and IGF-1R, respectively. In addition to ligand stimulation, cells were also treated with ligand blocking antibodies: MM-151, an oligoclonal therapeutic composed of three monoclonal antibodies targeting EGFR;23 Seribantumab (MM-121), a monoclonal antibody targeting ErbB3;16 MM-131, a bispecific antibody co-targeting Met and EpCAM;24 and Istiratumab (MM-141), a bispecific antibody co-targeting IGF-1R and ErbB3.25 Fig. Fig.1a1a illustrates the RTKs, their corresponding ligands and the mechanism of action of the ligand blocking antibodies. Proliferation was quantified in a 3D spheroid formation assay at the 3-day time point (Fig. (Fig.1b)1b) by measuring ATP content as surrogate for cell number (CellTiter-Glo® assay). Response was classified as positive if the signal at the 3-day time point was more than 20% above the respective control, plus being significant at a confidence level α=0.05 (measured in quadruplicates, Wilcoxon rank-sum test). Per this screen approximately 45% of cell lines responded to EGF, 55% of cell lines responded to HRG, 33% of cell lines responded to HGF and 7% of cell lines responded to IGF-1. The low response rate to IGF-1 in this proliferation screen may reflect the presence of IGF-1 in the low-serum medium and the modest absolute inhibition point to the importance of IGF-1 mediated signaling for survival rather than for proliferation.26 We and others observed a generally weaker MAPK activation via IGF-1R (see Fig. Fig.3c)3c) compared to the other growth factors in the screen.27,28 Further, the CellTiter-Glo® assay relies on metabolic function and hence can be limited as readout for IGF-1 stimulation.29
Figure Figure1b1b shows the response to treatment with ligand in combination with the respective blocking antibody compared to the ligand effect alone. Depending on the ligand treatment, 5–17% of cell lines were ligand non-responsive, but the antibodies inhibited basal proliferation, which is indicative of autocrine driven proliferation. Even though IGF-1 did not induce a proliferative response in most cell lines, MM-141 inhibited proliferation in about 19% of the cell lines indicating that IGF-1 might be present in low-serum medium.
Investigation of correlations between the ligand and antibody responses across all cell lines revealed a checkerboard pattern of significant positive correlations between EGF, HRG, and HGF as well as anti-correlations of those ligands and their respective antibody responses (Fig. (Fig.1c).1c). This suggests a general trend that cell lines are either responsive to multiple ligands and their respective antibodies (right hand side of Fig. Fig.1b),1b), or are generally non-responsive to any given ligand or antibody (left hand side of Fig. Fig.1b).1b). For IGF-1/IGF-1R, the only significant correlation was observed between Istiratumab treatment and Seribantumab treatment. This can be attributed to both antibodies (co-) targeting ErbB3 and, therefore, some cell lines respond to both Istiratumab and Seribantumab independent of an IGF-1 effect (see, e.g., KYSE-410 cell line in Fig. Fig.1b).1b). The general lack of correlation patterns for IGF-1/IGF-1R responses as were observed for the ErbB family and HGF/Met can be explained by the lack of IGF-1 induced proliferation in this screen.
In the following, we will focus on the question of how ligand dependence can be predicted. A necessary condition for response to any given ligand is the presence of its respective receptor. First, we used a univariate analysis (Fig. (Fig.1d)1d) and found that receptor surface levels measured by qFACS do not correlate significantly with the respective ligand response. Based on this data, a simple linear model cannot stratify responsiveness. Next, we investigated whether a multi-pathway signaling model featuring the complex receptor interactions as well as the cross-talk between the mitogen-activated protein (MAP) kinase and the phosphoinositide 3-kinase (PI3K) signaling pathways can be used to predict the phenotypic response. Specifically, signaling features like the area under curve (AUC), quasi steady-state and the signal amplitude of receptor homo- or heterodimers and downstream components were considered as inputs into a decision tree classification algorithm.
To construct a comprehensive signal transduction model that could be used to predict proliferation in response to growth factors for all 58 cell lines, we built on a previously published model of ErbB receptor signaling.16 We extended the computational model to include IGF1-R and Met (Fig. (Fig.2a)2a) as well as 12 homo and heterodimers for which biological evidence can be found.30–33 Our analysis considers EGFR, HER2, ErbB3, Met and IGF-1R homodimers as well as the heterodimers EGFR-HER2, EGFR-ErbB3, EGFR-Met, HER2-ErbB3, ErbB3-Met, IGF-1R-IGF-1R, EGFR-IGF-1R, and HER2-IGF-1R. The latter two were later removed from the computational model without impacting the model performance. Figure Figure2b2b depicts the structure of the model for the example of a signaling HER2-ErbB3 heterodimer. The complete model consists of 62 differential equations and replicates the model structure shown in Fig. Fig.2b2b for each of the considered ten homo- and heterodimers. In short, receptors bind ligand with published dissociation constants (K D). Bound receptors can form homo and heterodimers and subsequently undergo endocytosis. After internalization, the receptor dimers can get either dephosphorylated and recycled to the cell surface or they get degraded in the lysosomes.34,35 Downstream of the receptor, all homo- and heterodimers except the ErbB3-homodimer, which cannot trans-phosphorylate due to its lack of intrinsic kinase activity, can activate the MAP kinase cascade as well as the PI3K/AKT pathway. ERK and AKT phosphorylation converge in the phosphorylation of S6K1 and S6. Several known feedback mechanisms between the pathways 27 were implemented in the computational model. Mathematical details, executable code to simulate the model and instructions to replicate our findings are available in the supplementary materials and on biomodels.org.
The computational model is constructed with the aim to capture the signaling dynamics of key components of the signaling pathway including receptor homo- and heterodimerization. It is not intended to be a complete compendium of all the known molecular interactions.33,36,37 Size and complexity of the computational model were chosen to reflect the available experimental data, and to facilitate efficient computation. This is particularly important during model parameter calibration, which uses parameter estimation algorithms to match the available experimental data as closely as possible (see Methods section for details). Mutations were not implemented in the computational signaling model as they appear to increase the signaling baseline but not necessarily the signaling dynamics.38 However, the mutation status for each cell line was used in the machine learning classification.
For model calibration, phosphoproteomic time course data from protein microarrays39 for the receptor phosphorylation as well as for phospho-MEK, phospho-ERK, phospho-AKT, and phospho-S6 across all seven cancer cell lines (H322M, BxPc-3, A431, BT-20, ACHN, ADRr, and IGROV-1) were used. Only the two cell lines H322M and IGFROV-1 were included in the cell line proliferation screen in Fig. Fig.1.1. These seven cancer cell lines represent different cancer indications (lung adenocarcinoma, pancreatic, epidermoid, breast and ovarian cancer) and were selected based on the molecular diversity with respect to the mutation status and differences in receptor expression. A key challenge for building computational models that can describe and predict signaling dynamics of different cell lines is to limit the number of model parameters that are specific to one cell line.40 In this case, it was possible to restrict all kinetic rate constants to the same value and to adjust only the receptor expression for individual cell lines. Due to the analytically calculated basal activation levels of all homo- and heterodimers as well as of the downstream components, which were derived from steady-state constraints,41 the receptor expression impact the signaling response throughout the model. Therefore, the individual receptor expression of each cell line enables distinct model responses upon ligand stimulation. The receptor expression was measured using quantitative flow-cytometry (qFACS) in combination with RNA sequencing data (see Methods section). The model can accurately describe the time course data of seven training cell lines, with 85.7% of the data points within two standard deviations of the model uncertainty (see Fig. Fig.33 for a selection of the data and Suppl. Figs. 11–40 for a comprehensive comparison of model simulations and experimental data).
Based on the trained model, predictions were generated for two independent validation cell lines (BT-474 M3, MDA-MB-231) and compared to the experimental data. The goodness of the predictions for the validation cell lines was equivalent to the goodness of fit of the training cell lines (Fig. (Fig.3b,3b, Supplementary Figs. 35–40 for model fits to the available data). These simulation results validate that receptor expression is sufficient to predict signaling features of independent cell lines that were not used for model training. In addition, we generated model predictions that were based on random receptor surface levels, by taking non-matching values from randomly selected cell lines used in the cell viability assay (see Table Table3).3). The decline in goodness of fit was on average 30% and statistically significant (p=8.5 * 10−9, see Supplementary Fig. 2). These results illustrate the importance of receptor expression and their ratios to capture the distinct signaling features observed for each cell line.42
To further validate the presented model structure with its multiple receptor heterodimers, a simplified model lacking any heterodimerization capabilities was trained to the experimental data. In this setting, all receptors could signal downstream through homodimerization, disregarding the non-functional kinase unit of the ErbB3 receptor. Even with 39 parameters less, the reduced model had a goodness of fit impediment with associated p-value<1.e-15 in the corresponding likelihood-ratio test, showing the significant improvement of the computational model by including receptor heterodimerization. Besides, concordance of the basal receptor levels obtained via analytic steady state equations, reflecting the proposed receptor trafficking, was assessed through extensive measurements of basal total and phosphorylation levels in 39 breast cancer cell lines.28 A good correlation, especially for the ErbB receptor family, was found (see Supplementary Fig. 3), confirming the calibrated model parameters constituting cell-dependent steady states.
To further test the robustness and applicability of the model to ligands not included in the original training-set, we compared the predicted receptor activation patterns in response to different ligands of the EGFR-ligand family, such as Betacellulin (BTC). To this end, previously published16 time-resolved data of the ADRr cell line for EGF and BTC with ligand concentration range between 0.1nM and 10nM was reanalyzed with the current model (Suppl. Fig. 4a). Differences in the ligand binding affinities of each ligand to the EGF receptor as well as different homo- as well as heterodimerization kinetics were sufficient to describe the experimental data (Supplementary Figs. 4b, c). BTC induces a stronger EGFR homodimerization compared to the stronger EGFR-HER2 heterodimerization induced by EGF (Supplementary Fig. 4d). These differences in EGFR homo and heterodimerization with HER2 were previously described.43,44
Further insights into growth factor signaling and signal processing by the cancer cells can be gained by analyzing the computational model and why it can capture the distinct signaling dynamics across cell lines. This analysis revealed the importance of different homo and heterodimers in encoding information as a function of the ligand(s) present. The largest effect of heterodimerization on signal output is seen within the ErbB family. The interplay between the receptors explains the slower, more sustained receptor activation in response to EGF in the BxPc-3 cells, which are characterized by a high ratio of EGFR to other receptors (Fig. (Fig.3a).3a). In contrast to the BxPc-3 cells, the IGROV-1 cells are characterized by low EGFR expression compared to other receptors leading to the observed transient and early activation of EGFR and HER2. For the ACHN cell line, we generated signaling data in response to EGF, HGF as well as to the combination of EGF and HGF. The computational model reveals sophisticated feedback regulation between the MAPK and PI3K pathways, e.g., reduced AKT activation comparing EGF and HGF co-stimulation to HGF only (Fig. (Fig.3b).3b). In Fig. Fig.3c3c another example is depicted: when EGF and HRG are present, EGFR and ErbB3 compete for HER2. The ligand combination results in reduced phopsho-ErbB3 levels due to a dominant binding of HER2 to EGFR in the presence of EGF (Fig. (Fig.3c).3c). We argue that mechanistic understanding of changes in receptor stoichiometry based on individual ligands or ligand mixtures can cause non-obvious signaling responses and is required to understand the ultimate phenotypic response.
IGF-1 is distinct from the other growth factors in our screen. IGF-1 displayed a much weaker ability to induce proliferation and similarly the effect of IGF-1 signaling in co-stimulation experiments is distinct to the HRG/EGF or HGF/EGF co-stimulation experiments. The time-course data for co-stimulation of IGF-1 with either EGF or HRG did not show deviations from the respective stimulation with IGF-1 alone (see Supplementary Fig. 14). Consequently, the model parameters referring to the heterodimerization of IGF-1R (see heterodimers involving IGF-1R in Fig. Fig.2a)2a) with other receptors could be set to zero without a significant decline in the goodness of fit.
The calibrated mechanistic signaling model can be used to simulate signaling features for the stimulation with EGF, HRG, HGF or IGF-1 for all cell lines of the cell viability screen only using their receptor expression as inputs. To connect signaling features derived from the computational model to the phenotypic response observed in the cell proliferation screen, we applied a machine learning approach. Based on different sets of input features that were selected based on their prediction ability (see below), we trained bootstrap-aggregating (bagged) decision trees (BDTs). BDTs are highly efficient for multivariate analysis and allow for a comprehensive interpretation of the chosen features.45,46 Therein, a multitude of trees are trained, with each single tree aiming at discriminating growing from non-growing cells based on the provided feature space. At each node, a tree divides the data through selected features that yield the best improvement in signal-to-noise. Combining the ensemble of trees, high-dimensional and non-linear regions in feature space prone for cell growth are obtained and can be used for prediction. More details can be found in the supplementary materials and in Supplementary Fig. 5.
To investigate the contribution of dynamic signaling features derived from the computational model on the predictive performance of the machine learning approach, we considered two different sets of input features (additional sets are reported in Supplementary Fig. 6). The first feature set contained the receptor expression, KRAS and PI3K mutation status and ligand treatment as binary input (Fig. (Fig.4b4b).47,48 The second feature set did not contain the receptor expression explicitly but only signaling features derived from the computational model based on the receptor expression and the ligand stimulation (Fig. (Fig.4c)4c) as well as the mutation status (see Table Table3).3). The cell line specific signaling features consist of the AUC of all phosphorylated receptor homo- and heterodimers in addition to AKT, ERK, and S6 phosphorylation. Inclusion of the fold-change of these features as well as their quasi steady-state levels were tested but did not yield substantial benefits on top of the information given by the area under curve (see Supplementary Fig. 6). Figure Figure44 illustrates both multivariate prediction strategies as well as the univariate analysis shown in Fig. Fig.1d1d.
To evaluate the accuracy of BDT predictions, the cell lines were randomly split 500 times into training and testing sets. For each ligand, BDT training was performed on the training data for all available ligands, while efficiency was calculated on the testing cell lines for the chosen ligand only. By leaving out whole cell lines as opposed to a fraction of the total data, possible bias due to correlated responses to different ligands in the same cell line is avoided. We monitored the fraction of true predictions as a metric for the prediction accuracy. Both feature sets resulted in a better prediction of proliferation compared to random data, which results in 50% true predictions and serves as control (Fig. (Fig.5a).5a). Exceptions were the predictions for IGF-1 and HGF stimulation using the receptor expression only, where the performance drops insignificantly below the control. Training on features derived from the computational model improved the prediction of cell proliferation significantly compared to control (p-value of <1.e-2, see Fig. Fig.5b)5b) for the combination of all ligands except IGF-1, while BDT predictions based on receptor expression alone did not result in a statistically significant improvement (p-value of 0.15, see Fig. Fig.5b).5b). The respective distributions for single ligand induced proliferation predictions can be found in Suppl. Fig. 7. The BDT training was robust with respect to the relative amount of training and testing data and to the significance threshold upon which a cell is labeled as proliferating (Supplementary Fig. 8). For IGF-1, the low number of responders resulted in a low correlation within the events and a statistical bias of their relative amount in training or testing data. These circumstances rendered robust prediction impossible and the IGF-1 data set was excluded, e.g. in the combination of all ligands (see Fig. Fig.5a5a).
One of the advantages of mechanistic computational models is that it is possible to gain insights into cellular signal processing. Thus, the importance of different model features during BDT training can be traced back to develop hypotheses about what ultimately drives proliferation. The training features are ranked by their impact on data classification, which is measured by the average gain in signal-to-noise ratio over all trees. As illustrated in Table Table1,1, the most important features rely on the homo- and heterodimerization stoichiometry of the ErbB receptor family as well as on the downstream signaling. A more detailed overview is given in Fig. Fig.5c,5c, which illustrates the data from the proliferation screen and the model-derived features that proved important during training of the BDT (see Table Table1).1). This coincides with prior biological knowledge that ErbB receptors induce proliferation.13 It can be observed that EGF, HGF or HRG stimulation induce very specific homo and heterodimerization patterns of EGFR and EGFR/HER2 respectively. Moreover, phospho-S6 is an important feature to predict proliferation in the presence of HRG and HGF, both mainly activating the PI3K pathway. Its importance might be a result of the crosstalk between the MAPK and PI3K pathways and the convergence of both pathways in S6 phosphorylation.49 Moreover, PI3K mutation status and EGFR heterodimerization patterns help to identify clusters of EGF dependent cell lines. RAS mutation status and differences in activation of downstream targets further predict proliferation after HRG and HGF stimulation. Independent of the KRAS mutation status, HRG induces increased levels of AKT phosphorylation while inducing the same amount of S6 phosphorylation as HGF.
In the previous section, we showed that a mechanistic computational model in combination with decision tree classification can predict in vitro proliferation. Next, we applied this novel approach to patient data. Using the data from the TCGA Research Network (http://cancergenome.nih.gov/), we use our model to predict if an individual patient tumor would show a proliferation response if stimulated by the ligand of interest. The patient data set includes 2909 samples from patients with breast, colorectal, lung, and ovarian cancer. The input to our model is receptor RNA expression measured by RNA sequencing for each tumor sample. The measured RNA expressions between our in vitro cancer cell line data and the data from patient tumors were on different scales. Therefore, we normalized the expressions to their respective means. Subsequently, the expression of EGFR, Her2, ErbB3, Met, and IGF-1R were used to perform model simulations and to extract the signaling features needed to predict ligand-dependent tumors using the previously described BDT algorithm. The number of ligand-dependent tumors differed within indications and ligand (EGF, HGF or HRG). The number of predicted ligand responsive tumors was highest for HRG followed by EGF and lowest for HGF (Fig. (Fig.6a).6a). Lung and colorectal cancer seem to be most responsive to EGF, which is congruent with the high prevalence of EGFR mutations and overexpression in these indications.50–52 In Fig. Fig.1b1b we observed that ligand induced proliferation is correlated with treatment response to an antibody targeting the respective receptor. Therefore, the approvals of EGFR inhibitors in non-small cell lung cancer (NSCLC) and CRC confirm the predicted dependence on EGFR signaling.53,54 In contrast, the low dependence of NSCLC on HGF signaling might explain the failure of Onartuzumab (MetMAb), a Met blocking antibody in a Phase 3 study in NSCLC.55 Similarly, EGFR inhibitors have not yet proven to result in clinical benefit in breast cancer,56 which is also in agreement with the predicted low EGF dependence of breast cancer. The predicted high responsiveness of breast cancer and lung cancer to HRG seems to agree with the retrospective analysis of two clinical studies with Seribantumab and the finding that HRG expression appears to be predictive of patients responding to therapy.57 Unfortunately, the precision of the predictions cannot be assessed as no outcome data to treatment with ligand blocking antibodies is available for the TCGA data set. The predicted ligand-dependence only considers the molecular makeup of individual tumors. However, a tumor that is predicted to be ligand-dependent would respond only if the respective ligands were also present in sufficient amounts in the tumor microenvironment. The local concentration of ligands, however, cannot be inferred from our analysis and it is difficult to match to the in vitro data that was used for model training. This impacts the predicted number of ligand-dependent tumors in this data set.
However, apart from the relative number of ligand-dependent tumors, we observed a significant correlation between ligand expression and the predicted response to ligands (Fig. (Fig.6b,6b, additional data in Supplementary Fig. 9). The predicted ligand-dependent tumor samples from patients with breast and colorectal cancer display statistically significant higher (t-test) amounts of the corresponding ligand compared to the predicted ligand independent tumors, if we compare the mean expression. This suggests that tumors that express ligands evolved to be sensitive to ligands, or vice versa. In our opinion, this is an indirect proof that the model predictions can be applied to data from patient tumors and that they could be clinically relevant.
The research to understand and therapeutically battle various cancer types has made significant progress over the last two decades, decreasing the overall mortality by roughly 2% per year since 2001.58 Targeted therapy and combinations thereof have become important areas of drug development with rising FDA approval rates in the last years.59,60 However, by studying cellular fate across multiple cell lines and indications, we and others have learned that complex interactions between receptors as well as positive and negative feedback regulation between signaling pathways can diminish drug efficacy. To obtain a deeper understanding of cell response to exogenous stimuli, phenotypic responses need to be studied in the context of multiple signaling pathways as well as mutation status. In this work, we developed a computational model describing multiple signaling pathways and show that a BDT algorithm using simulated signaling features can accurately predict ligand-dependent proliferation in vitro.
The signaling model incorporates the ErbB receptor family as well as the Met and IGF-1R receptors. Parameters of the model were estimated based on a variety of time-resolved data from seven different cell lines including a wide range of ligand concentrations with comprehensive single ligand and co-stimulations. The cell lines cover a broad range of ratios of receptor expression. The ligand concentrations used in the proliferation screen were in the range of concentrations used for the signaling experiments. While retaining a good fit to the experimental data, we could keep all kinetic parameters in the signaling model constant and just vary the receptor expression to describe the experimental signaling data for each cell line. It is important to note that this is not a direct proof that the reaction rate constants are identical between the cell lines. However, our model with all rate constants set to the same values, is in line with the receptor phosphorylation and downstream signaling data. We do not argue that the model presented here is correct in all of its aspect, but we could show that the signaling dynamics predicted by the model are useful for predicting the cellular responses to ligand stimulation and that this presents an approach that should be explored further by incorporating clinical patient response data. The computational model not only describes the data for the seven training cell lines but also predicts the signaling responses of two additional, independent validation cell lines. We showed that the goodness of fit is dependent on the absolute receptor expression and the formation of different homo- and heterodimers. The observed variability in model response was achieved by differences in internalization, degradation and recycling rates for different receptor homo- and heterodimers. Saturation of different downstream model components is sensitive to the receptor expression. The analytically solved steady states for all cell lines were important to implement the complex receptor dimerization properties in the signaling model. Thus, ligand stimulation results in a complex re-distribution of receptor homo- or heterodimerization and facilitated distinct downstream activation patterns. The calculated steady states, especially for the ErbB receptor family, were found in concordance with measured basal total and phosphorylation levels of 39 breast cancer cell lines utilized from.28 The fact that the developed computational model can accurately describe the signaling responses across multiple cell lines enables the prediction of signaling dynamics for cell lines of the proliferation screen that were not used to construct the computational signaling model.
To predict proliferation in response to ligand stimulation, we linked simulated signaling features to the data of the cell proliferation screen using a supervised machine learning approach. Tree-based classification algorithms are widely used for machine learning61–63 and carve out regions in feature space that best distinguish between different data classifications, here proliferation vs. stasis upon ligand stimulation. BDTs were trained on the in vitro proliferation screen across 58 cell lines. The algorithm was trained with either receptor expression and ligand stimulation conditions as Boolean columns or with signaling features extracted from the computational model. The signaling features included the integrated area under curve of all receptor complexes as well as the phosphorylation of AKT, ERK, and S6. Both feature sets lead to a better prediction of proliferation compared to control, where response was predicted based on random data. However, the signaling features allowed for a more robust and statistically improved prediction of proliferation. In addition, the computational model allows us to gain insights into the underlying processes driving ligand-dependent proliferation. In all cases homo or heterodimerization of the ErbB receptors was important. In the case of EGF, the PI3K mutation status mattered in addition. For HRG and HGF possible RAS mutations together with AKT and S6 phosphorylation were important features to predict cell proliferation. However, the importance of features in the tree-based approach is very sensitive to the utilized data, and additional measurements are needed to infer the role of MAPK and PI3K signaling in inducing growth.
Simulated signaling features are advantageous over using receptor expression directly as input features for two reasons: First, the dynamic range of receptor activation as well as of the downstream components is described quantitatively and renders the model outputs more robust to receptor expressions, which span multiple orders of magnitude. Second, the interplay between receptors and the included feedback mechanisms adds a source of information on top of the receptor expression and ligand information alone, resulting in a non-linear input transformation that improves the detection of regions governing proliferation.
To demonstrate the applicability of this novel approach to patient samples, data from 2909 patients with breast, colorectal, lung or ovarian cancer were analyzed. For these samples, model simulations were conducted to extract signaling features required for BDT prediction of ligand-dependence. Interestingly, we observed a significant correlation between the measured ligand expression and the predicted ligand-dependence for breast and colorectal cancer. This may be a consequence of evolutionary adaption of these tumor cells due to the growth advantage from ligand-mediated signaling. Therefore, the presence of ligands in the tumor micro-environment may be a favorable biomarker for RTK-directed drug treatment. This rationale together with possible mutations, which contributed substantially in the BDT training, are currently being explored in the clinic.57
However, both the prediction of proliferation and the mechanistic computational model have limitations. For one, machine learning in feature-space regions that are not covered by many cell lines is not efficient as illustrated by IGF-1 in the proliferation screen data-set. The mechanistic model is limited in its ability to fully reproduce data in cases of either receptor overexpression (see Supplementary Fig. 19), which probably transfers to the presence of activating mutations, as well as in cell lines harboring PI3K or RAS mutations (see Supplementary Figs. 22, 23, and 37 to 40). We also encountered computational limitations since the complexity of the mechanistic computational model is on the verge of what is currently computationally feasible. This said, more emphasis on model selection, e.g. profile likelihood-based model reduction64 with additional prior knowledge may allow us to better bridge between quantitative time-resolved data and large-scale genomic and phenotypic data. To understand ligand mixtures and pathway redundancy a greater variety of single ligands, e.g., FGF and PDGF, or stimulations with ligand mixtures might aid to more accurately determine model parameters for receptor dimerization, trafficking and downstream activation in the future. Further, or model is currently limited to MAPK and PI3K pathways. Future work, should expand on this approach by including other signaling pathways, e.g., data from the JAK-STAT signaling pathway. An expanded model may improve the accuracy of the predictions as well as additional perturbations like specific gene knockdown etc. would help to improve the signaling model.
Ligand-dependence or addiction to growth factors (ligands) might prove to be far more prevalent than oncogene addiction given the high ligand prevalence in solid tumors and potentially as much more complex since multiple ligands are expressed in the tumor microenvironment. To successfully treat patients with ligand-dependent tumors with targeted inhibitors (small molecules or monoclonal antibodies) or rational combinations of targeted inhibitors, a better understanding of ligand-dependence is crucial, especially given the redundancy of signaling pathways within a tumor cell and tumor heterogeneity. We argue that the mechanistic understanding of changes in receptor stoichiometry based on individual ligands or ligand mixtures result in non-obvious signaling responses that are relevant to the ultimate phenotypic response. The work presented here demonstrates that for targeted therapies to be successful in the clinic the ligand hierarchy as well as co-dependence need to be understood and most likely require the measurement of multiple ligands and respective receptor expression in tumor biopsies. New methodologies like single cell RNAseq65 may allow us in the future to characterize the clonal composition of tumors and to determine which cellular fraction is ligand-dependent and which drug combination is best suited to eliminate all ligand-dependent tumor cells.
The presented novel approach of using BDTs in conjunction with simulated signaling features is the beginning of how complex mechanistic models and large data sets can be combined to understand cell-specific complexity but also heterogeneous tumors better. We demonstrated that mechanistic computational models of signaling pathways can help bridge between large scale in vitro observations and clinical hypotheses. In the future, selected in vivo studies should be used to validate rational combination regimens. Previous efforts predicting drug sensitivity based on large and diverse data sets found that gene expression data proved most valuable, together with exploiting non-linear relationships and addition of prior knowledge of biological pathways.18,66,67 Yet, significant improvement in predictions proved to be challenging across multiple approaches and data sets. With the presented approach, mechanistic knowledge can be easily combined with known datasets from RNA sequencing, copy-number alterations and mutation information to improve the prediction of patient-individual drug response and unravel the interplay between complex signaling and cellular fate.
All cell lines used in the viability screen were purchased from American Type Culture Collection. In brief, cells were cultured in RPMI 1640 medium (Life Technologies) supplemented with 10% fetal bovine serum (FBS, Life Technologies) and 1% penicillin-streptomycin (pen/strep, Life Technologies) at 37°C and 5% CO2. Recombinant human heregulin1−β1 (NRG1b) EGF domain and HGF were obtained from PeproTech, and recombinant human IGF-1 and EGF were obtained from RD Systems. Seribantumab, MM-131, Istiratumab and MM-151 was manufactured in-house by the Merrimack Pharmaceuticals Protein Engineering Department and stored at 4°C. CellTiter-Glo® was obtained from Promega and reconstituted fresh for each experiment.
To measure cell viability in a three-dimensional spheroid culture, cells were seeded into 384-well low- binding multi-spheroid culture plates (Scivax, USA) in the relevant growth medium supplemented with 4% FBS and 1% pen/strep at a density of 1500 cells/well. To allow for spheroid formation, plates were incubated for 24h, after which cells were treated with ligands (5nM HRG1, 5nM EGF, 50nM IGF-1, 1nM HGF) and/or inhibitors (1μM of Seribantumab, MM-131 and MM-151 and 0.5μM of Istiratumab) in 4% FBS containing medium. Following 72h of incubation, cell viability was determined by incubation with CellTiterGlo® reagent for 10min and measuring well luminescence on an Envision (Perkin Elmer) plate reader.
The available data consists of three different data sets, which include time-resolved concentration measurements of activated and total receptors as well as of various phosphorylated downstream targets. The largest data set comprises nine cell lines with measurements of the EGFR, HER2 and ErbB3 receptors of the ErbB family and the IGF1-receptor, together with the downstream targets ERK, AKT, S6K1 and S6. Four different ligand concentrations of EGF, HRG, and IGF-1 ranging from 0.156 to 10nM are used and 12 measurement time points up to 240min are taken. In addition, co-stimulations of the respective ligands are available for two of the nine cell lines. Out of the nine cell lines, six cell lines are used for model calibration, the remaining three to validate the model. Cell lines used for calibration include H322M (non-small cell lung cancer), BxPc-3 (pancreatic cancer), A431 (epidermoid cancer), BT-20 (breast cancer), ADRr (ovarian cancer) and IGROV-1 (ovarian cancer). BT474 (breast cancer), MDA-MB-231 (breast cancer) and ACHN (renal cancer) are utilized to validate the model.
Measurements after either EGF or BTC stimulation are available for one of the calibration cell lines, ADRr, with ligand concentrations between 0.11 and 9.26nM, and 12 measurement time points up to 240min. These include phosphorylation of EGFR, HER2, ErbB3, ERK, and AKT. Apart from that, measurements with HGF and EGF as well as their co-stimulation is available for ACHN, which is used for validation with respect to HRG and IGF-1, spanning the same concentrations and measurements up to 120min. Therein, phosphorylated EGFR and Met phosphorylation as well as phospho-ERK and phospho-AKT are measured. The receptor concentrations in all experiments are measured by ELISA whereas the downstream components are measured by lysate microarray.
Cells were trypsinized, washed, and stained using fluorescently labeled antibodies, see list below. Antibodies were labeled as previously described (16). Receptor numbers were determined by assessing the antibody-binding capacity of the fluorescently labeled antibody via quantitative fluorescence-activated cell sorting. Antibody-binding capacity was determined using Simply CellularQuantumBeads (BangsLabs, Fishers, IN), per the manufacturer’s instructions. List of antibodies and target receptors used in qFACS method: Erbitux (EGFR), Trastuzumab (HER2), Anti-human HER3 Ab generated by Merrimack, A12 IgG (anti-human IGF1-R), Anti-human Met, Mouse Anti-Human EpCAM-APC (BD Biosciences, Cat# EBA-1).
Matching high density lysate matrices from nine cells lines (BxPc-3, H322M, ACHN, IGROV-1, A-431, BT-20, ADRr, BT-474-M3, and MDA-MB-231) were generated for both the ELISA and lysate microarray studies. Lysate matrices were developed by treating each cell line with EGF, HRG1b1, IGF1, individually at four different doses (10nM, 2.5nM, 0.625nM, and 0.156nM), as well as in all two-way combinations of these three ligands at a single 2.5nM dose of each ligand (a total of 15 conditions) in 96-well plates (Corning). Antibodies directed at pAKT (S473), pMEK (S217/S221), pERK (T202/Y204), p-S6K1 (T389), and pS6 (S235/S236) were obtained from Cell Signaling Technology. Cells were cultured using standard tissue culture techniques in RPMI supplemented with 10% FBS and penicillin/streptomycin. Cells were counted using a hemocytometer and seeded at 10,000 cells/well into 40x 96-well tissue culture plates (Corning). After 24–48h, once cells had reached approximately 50% confluency, medium was aspirated and replaced with low- in medium (RPMI supplemented with 0.5% FBS and penicillin/streptomycin) (Gibco). 16–24h later, 2x ligand solutions prepared in low-serum medium were added simultaneously to all 96 wells of each plate. Plates were returned to the incubator for the prescribed incubation times, then placed on ice to stop ligand stimulation. Medium was aspirated from all wells, cells were washed once with ice-cold PBS, and then lysed in 30μL/well of lysis buffer. Twelve time points as well as an untreated time zero were collected for all conditions (0, 2, 4, 6, 8, 10, 15, 30, 60, 90, 120, and 240min). Lysates were then collected and frozen. For ELISA lysates, M-PER buffer (ThermoFisher) was supplemented with protease and phosphatase inhibitors (Roche) and NaCl to a final concentration of 150mM. For lysate microarrays, lysis buffer was prepared as previously described.68
Total and phospho ELISAs of EGFR, ErbB2, ErbB3, and IGF-1R were performed as previously described.69 Protein specific antibodies were used for capture in all cases, while an anti-phospho tyrosine antibody was used for detection in all phosphor assays (see Supplementary Table 1). Antibody screening for this study and well as lysate preparation and printing were performed as previously described.68,70 Arrays were printed with a 7-point dilution series for each lysate onto 16-pad slides (GraceBioLabs OnCyte Avid, Bend, OR) using the Aushon 2170 micro-arraying robot in 2x4 8-pin mode (Aushon Biosystems, Billerica, MA). Slides were washed in 100mM Tris-HCl pH 9.0 at RT for several days, followed by 3x 5min PBST, after which slides were spun dry. ProPlate 16-well slide modules (GraceBioLabs) were then attached. Arrays were blocked in Odyssey Blocking Buffer (LiCor, Lincoln, NE) at 4°C for 1h, after which blocking solution was replaced with 1:500 primary antibodies (all rabbit, see list below) mixed with 1:1000 anti-beta-actin antibody (mouse) (Sigma A1978) for 24h at 4°C with agitation. Slides were then washed 3x 5min with PBST, then incubated in secondary antibodies (1:1000 anti-rabbit-800 and 1:1000 anti-mouse-680 (LiCor) in 5% BSA/PBST) at 4°C for 5h. Arrays were washed briefly in PBST, ProPlate modules were removed, and whole-slides were washed 3x 5min in PBST. Slides were then spun dry at room temperature. Slides were scanned on the LiCor Odyssey scanner at 21μm resolution and under the “highest” quality setting in both the 700 and 800nm channels. Spot intensities were extracted using the LiCor ImageStudio software with manual spot alignment.
Mechanistic models based on ordinary differential equations (ODEs) are frequently used for the description of biochemical reaction networks. They are composed of kinetic rate equations and every component x of the model has a biological counterpart. The time evolution x(t) of the model concentrations is obtained by integration of the corresponding system of ODEs
depending on initial and kinetic rate parameters comprised in θ d. These are linked to measured concentrations of the involved constituents y(t) by an observational function
with the assumption of Gaussian errors ϵ~N(0, σ) that is often achieved via log transformation. In addition, the observation function includes e.g. scaling and offset parameters, summarized in θ o. Both observational and dynamic parameters are comprised in θ. To compare the model response to measured data at time points t i, the scaled log-likelihood is calculated via
Within the maximum likelihood framework, the optimized parameter set is estimated through minimization of χ 2 (θ).
Since analytical solutions of non-linear ODE systems are in general not available, a numerical integration must be performed. In this work, the dynamical system and its sensitivities were integrated by the CVODES integrator of the SUNDIALS suite.71 Therein, an implicit BDF integration method72 with attached KLU sparse solver was chosen.73 The inner derivatives of the likelihood needed in gradient-based parameter estimation were computed via forward sensitivities supplied to the integration algorithm.74 Numerical optimization was conducted using a trust-region based, large scale nonlinear optimization algorithm implemented in the MATLAB function LSQNONLIN.75 For the mathematical modeling and visualization, the open-source and freely available d2d framework,76 based on MATLAB, was used.
We established a relationship between receptor levels on the cell surface measured by qFACS (see Methods section) and receptor mRNA expression from the Cancer Cell Line Encyclopedia (CCLE) for 124 cancer cell lines.77 A good correlation between mRNA and protein expression was previously shown in an independent study.78 By fitting a linear model (see Suppl. Fig. 10) we could calculate receptor surface levels or receptor mRNA expression also for the cell lines where data was missing (see Table Table2,2, Table Table3,3, and Supplementary Fig. 41).
The signaling model and code for machine learning analysis from this publication have been deposited to BioModels.org with the identifier MODEL1708210000. In addition, the computational model including all proteomic data, the phenotypic data from the in vitro viability screen together with the RNA sequencing data obtained from the CCLE 77 and the TCGA Research Network (http://cancergenome.nih.gov/), respectively, have been deposited within freely available modeling toolbox Data2Dynamics (http://data2dynamics.org; repository folder Examples/Hass_npjSysBio2017). It includes MATLAB code and documented script files to readily perform all analysis steps outlined in this publication. The code folder is also available in the Figshare repository (10.6084/m9.figshare.5331544). In addition, the mechanistic signaling model was also implemented in the open-source R package dMod 71 (https://github.com/dkaschek/dMod; Figshare repository: 10.6084/m9.figshare.5336338). The package contains the experimental data for all calibration cell lines and allows to simulate model trajectories.
We thank Tim Heinemann, Jeffrey Kearns, Sergio Iadevaia, Yasmin Hashambhoy-Ramsay, and Tim Maiwald for their constructive feedback and proof reading the manuscript, and Daniel Kaschek for implementing the model in R.
HH build the mechanistic signaling model, built BDT and TCGA predictions, performed the computational analysis and wrote manuscript. KM performed the cell viability screen and wrote the manuscript. MS and JA performed the ELISA and lysate microarray assays. VP performed the qFACS assay. SW helped built the mechanistic signaling model. JT and JS revised the manuscript. BS and GM helped plan the study and wrote the manuscript. AR planned the study, supervised computational work and wrote the manuscript.
The authors declare no competing financial interests.
Electronic supplementary material
Supplementary information accompanies the paper on the npj Systems Biology and Applications website (10.1038/s41540-017-0030-3).
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.