Predicting ligand-dependent tumors from multi-dimensional signaling features

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.


INTRODUCTION
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][4][5] Identifying these relatively rare patients via predictive diagnostic tests relying on genomic biomarkers has created Precision Medicine. [6][7][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][10][11][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][14][15][16][17] Large collections of genetic and genomic data led to efforts to disentangle the complex mechanisms using machine-learning algorithms. [18][19][20][21] It was previously shown that simulated patientspecific 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 liganddependent 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 coevolution of ligand production and the ability of the tumor cells to respond to stimulation.

RESULTS
In vitro proliferation screen 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. 1a 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 Results of the proliferation screen across 58 cell lines. Dots mark a significant increase in ligand induced proliferation or decrease in the presence of ligand plus antibody. The ligand effect is normalized to the medium control, whereas the antibody plus ligand effect is relative to ligand alone. The two cell lines marked with an arrow, as well as five additional cell lines that were not included in the proliferation screen, were used to train the computational model to signaling data. c Correlation pattern of ligand and antibody effects across all cell lines. d Linear correlation of receptor expression to ligand induced proliferation. The proliferation in response to ligand (y-axis) is displayed as log10-fold change with respect to day 0. The receptor surface levels (x-axis) are absolute measurements of receptors/cell by qFACS on a log10-scale point (Fig. 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. 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 1b 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. 1c). This suggests a general trend that cell lines are either responsive to multiple ligands and their respective antibodies (right hand side of Fig. 1b), or are generally non-responsive to any given ligand or antibody (left hand side of Fig. 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. 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. 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.
Multi-pathway computational model 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. 2a) as well as 12 homo and heterodimers for which biological evidence can be found. [30][31][32][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 2b 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. 2b 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 transphosphorylate 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 microarrays 39 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. 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. 3 for a selection of the data and Suppl.  [35][36][37][38][39][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 nonmatching values from randomly selected cell lines used in the cell viability assay (see Table 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 published 16 time-resolved data of the ADRr cell line for EGF and BTC with ligand concentration range between 0.1 nM and 10 nM 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 Importance of receptor homo and heterodimers 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.  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. 3b). In Fig. 3c 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. 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.  2a) with other receptors could be set to zero without a significant decline in the goodness of fit.
RNAseq and signaling features are predictive of phenotype 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 Fig. 4 Strategies for predicting ligand-induced phenotypic response. Based on the receptor expression of individual cancer cell lines, either a univariate or multivariate approach can be used to predict the phenotypic response to ligand stimulation. a Univariate approaches relate the respective receptor expression to the observed ligand induced proliferation for each of the four ligands separately. b-c Multivariate approaches such as bagged decision trees (BDTs) relate high-dimensional feature sets to the observed phenotype. b In this case the feature set consists of the five receptor surface levels as well as information about the respective ligand stimulation and mutation status. c The calibrated and validated signaling model allows to simulate the expected signaling dynamics for each individual cell line based on its receptor expression and ligands present. Based on the mechanistic knowledge that the signaling model incorporates, it can expand the initial fivedimensional feature set to a 12-dimensional feature set. This expanded feature set, together with information about mutation status is now connected to the observed growth responses by a bagged decision tree 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, highdimensional 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. 4b). 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. 4c) as well as the mutation status (see Table 3). The cell line specific signaling features consist of the AUC of all phosphorylated receptor homoand 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 4 illustrates both multivariate prediction strategies as well as the univariate analysis shown in Fig. 1d.
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. 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. 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.  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. 5a).
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 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. 5c, which illustrates the data from the proliferation screen and the model-derived features that proved important during training of the BDT (see Table 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.
Application to patient data 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. 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][51][52] In Fig. 1b 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 liganddependence only considers the molecular makeup of individual tumors. However, a tumor that is predicted to be liganddependent 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 Predicting ligand-dependent tumors H Hass et al.
expression and the predicted response to ligands (Fig. 6b, additional data in Supplementary Fig. 9). The predicted liganddependent 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.  5 Prediction of ligand-induced proliferation using BDTs. a Ratio of true predictions after BDT training with simulated signaling features or receptor expression only, compared to random predictions in the presence of EGF, HRG, IGF or HGF. b For 500 random splits of training and testing cell lines, the BDT outcome is compared to random growth assessment as histogram and cumulative density function, showing the significant improvement due to mechanistic modeling. c Data of in-vitro cell viability screen showing proliferation response (green) or no significant response (red) in different 2D representations of the feature space

DISCUSSION
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 costimulations. 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  Fig. 6 Predicting ligand dependent tumors from the TCGA data set. a Predicted percentage of tumors that would response to ligand exposure. The predictions were obtained by using the receptor RNA expression measured in breast, colorectal, lung, and ovarian cancers as inputs to our model. b The measured RNA expression of the ligands in predicted responders (red) vs. non-responders (green). The mean expression (black horizontal lines) and statistical significance of differences is indicated. The receptor mRNA expression is measured in transcripts per million and is displayed on a log-2 scale 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 homoor 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 learning [61][62][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 treebased 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 reduction 64 with additional prior knowledge may allow us to better bridge between quantitative time-resolved data and largescale 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 nonobvious 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 RNAseq 65 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 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, copynumber alterations and mutation information to improve the prediction of patient-individual drug response and unravel the interplay between complex signaling and cellular fate.  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 120 min. 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.

Cell lines and reagents
Quantification of receptor expression in cell lines using qFACS 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 antibodybinding 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), Antihuman Met, Mouse Anti-Human EpCAM-APC (BD Biosciences, Cat# EBA-1).

ELISA and Lysate microarray (reverse phase protein array) measurements
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 (10 nM, 2.5 nM, 0.625 nM, and 0.156 nM), as well as in all two-way combinations of these three ligands at a single 2.5 nM 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-48 h, 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-24 h 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 240 min). 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 150 mM. 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 100 mM Tris-HCl pH 9.0 at RT for several days, followed by 3x 5 min 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 1 h, after which blocking solution was replaced with 1:500 primary antibodies (all rabbit, see list below) mixed with 1:1000 antibeta-actin antibody (mouse) (Sigma A1978) for 24 h at 4°C with agitation. Slides were then washed 3x 5 min 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 5 h. Arrays were washed briefly in PBST, ProPlate modules were removed, and whole-slides were washed 3x 5 min 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 800 nm channels. Spot intensities were extracted using the LiCor ImageStudio software with manual spot alignment.

Mechanistic modeling
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.
Predicting ligand-dependent tumors H Hass et al.
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 method 72 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.

Calculation of receptor surface levels
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 2, Table 3, and Supplementary Fig. 41).

Data availability
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