Integrative Modeling Identifies Key Determinants of Inhibitor Sensitivity in Breast Cancer Cell Lines

Breast cancer is the most commonly diagnosed cancer in women and the second leading cause of cancer-related death (1). Decades of research have increased our knowledge of the molecular basis of this disease, whereas recent large-scale genomics studies have provided detailed information on mutations, copy number aberrations, and gene expression across different tumor samples, including breast cancer (2, 3). However, despite this increasing amount of knowledge, response rates for cancer treatment remain very low for many cancer subtypes, including breast cancer.

One of the challenges of cancer treatment is the genetic complexity of the disease, involving different oncogenic drivers or combinations thereof that allow the tumor to grow and proliferate. In some subtypes of breast cancer, the major oncogenic drivers are known, as is the case in HER2-amplified breast cancer, associated with the overexpression and aberrant activation of HER2 receptor signaling. Targeted treatment against HER2 indeed provides clinical benefit (4, 5). However, intrinsic resistance is frequently encountered, and acquired resistance often develops in tumors that are initially sensitive (6). A variety of drug resistance mechanisms to HER2-targeted therapies have been discovered in cell lines. For example, activation of the PI3K pathway resulting from mutations in PI3K (7), the loss of PTEN (8), or autocrine HGF signaling (9) have all been reported as mechanisms of drug resistance. However, due to the multitude of resistance mechanisms, which is further complicated by the cross-talk in downstream signaling, it is unclear to what extent each of these mechanisms is important for determining the sensitivity of a particular cell line or tumor. In other types of breast cancer, in particular triple-negative breast cancer, the regulatory signaling that drives the growth of the tumor is even less clear, although the PI3K/AKT pathway deregulation has been identified as a recurrent event (2).

Similar to patients, cell lines show a large degree of variability in drug response (10–12). Understanding the heterogeneous response in cell lines is an important starting point for understanding patient response. But despite many efforts, fully explaining drug response in vitro remains a challenge. Computational modeling can be used as a tool capable of untangling the complexity of drug sensitivity and resistance across different cell lines. There have been various approaches to computational modeling of drug sensitivity. These approaches can be broadly divided into two categories: approaches using linear or black box statistical models (11, 13, 14), and approaches using more detailed mechanistic computational models (15, 16). Although both approaches have recovered known drug sensitivity mechanisms and identified several novel associations, they have several limitations. For the black box statistical models, such as elastic net regression (10, 11), random forests (13), support vector machines (13), or a clustering-based method named ACME (14), the models are not sufficiently detailed to capture how molecular characteristics affect drug sensitivity. For example, the interactions between molecular aberrations are not explicitly modeled. This precludes finding all but the strongest associations, despite the very large number of cell lines that were profiled. In addition, available knowledge of signaling pathways is not employed, which could increase the statistical power to find molecular mechanism that associate with drug sensitivity. In the more detailed, mechanistic computational models (15, 16) such background knowledge is used; however, in these cases the number of cell lines studied has been limited, and it is thus unclear to what extent the particular mechanisms are important for explaining the variability across a larger set of cell lines. In addition, these mechanistic modeling studies used only a single data type, for example (phospho)protein expression, limiting the insight into the impact of other molecular aberrations present in the cell lines examined.

To address these limitations, we set out to combine detailed computational modeling of drug sensitivity mechanisms with extensive measurements of multiple data types derived from a breast cancer cell line panel. We developed a combined experimental/computational modeling approach, which can utilize background knowledge from the literature and integrate diverse types of data, including DNA sequencing, RNA sequencing, protein expression, and protein phosphorylation, with drug response data. We subsequently employed the computational model to analyze how the regulatory signaling in each cell line influences response to each drug. The computational model can also be used to identify cases where drug response cannot be explained fully by the existing knowledge. In one case, this lead us to identify and confirm the level of expression of eukaryotic translation initiation factor 4E-binding protein 1 (4E-BP1) as a determinant of response to mTOR inhibitors in breast cancer cell lines.

Together, we show the utility of employing integrative computational modeling to combine prior knowledge with measurements of multiple molecular data types to systematically test and extend our understanding of drug response to kinase inhibitors.

Detailed materials and methods are provided in the Supplementary Materials and Methods.

A panel of 30 breast cancer cell lines was assembled from various sources, the details and growth conditions of which are listed in Supplementary Table S1. Cell lines were authenticated by suppliers using short tandem repeat profiling. The SK-BR-7 cell line was obtained from an internal NKI cell bank and authenticated by STR profiling. We further confirmed the identity of the lines by comparing their mutation profiles (found by DNA sequencing) with those reported in COSMIC (17). Upon receipt, all lines were first expanded and early passage stocks frozen in liquid nitrogen. Lines were kept in culture for no more than 3 months, after which a new cell aliquot was obtained from frozen stock if needed. All cell lines were tested in-house and found to be negative for mycoplasma. Doubling times for each cell line were estimated by fitting exponential curves to confluence measurements obtained using an IncuCyte FLR/ZOOM instrument (Essen Bioscience). The resulting doubling times are provided in Supplementary Table S2.

Prior to carrying out drug response assays, cell line seeding densities were optimized. Cells seeded as for the cell doubling time experiments were assessed at the 96-hour endpoint for percentage confluence, and incubated with CellTiter-Blue (CTB; Promega) for a measure of metabolic activity. This was to ensure that cell lines did not exceed 90% confluence at assay endpoint and that the CTB signal at this density was not saturated. Seeding densities used for each cell line are listed in Supplementary Table S1.

For drug response assays, cells were seeded at the optimized density and 24 hours later treated with a 10-point 1:3 dilution series of a number of inhibitors using a Microlab STAR workstation fitted with 8 × 1,000 μL channels and 96-probe head (Hamilton): AZD8055, top dose 3 × 10⁻⁵ mol/L; BEZ235, 1 × 10⁻⁵ mol/L; GDC0941, 3 × 10⁻⁵ mol/L; MK2206, 3 × 10⁻⁵ mol/L; PD0325901, 3 × 10⁻⁵ mol/L; lapatinib, 3 × 10⁻⁵ mol/L; foretinib, 3 × 10⁻⁵ mol/L (all from Selleckchem). Each condition, including an untreated negative control and a phenyl arsine oxide (1 × 10⁻⁶ mol/L) treated positive control were set up in technical quadruplicate. Following 72-hour incubation, cells were stained with CTB (1:30 dilution) for 4 hours and the signal measured using an Envision spectrophotometer (Perkin Elmer). In the case of the validation experiments with HCC1806 and HCC1937 cell lines expressing 4E-BP1 or GFP constructs, cells were treated in a 9-point 1:3 dilution series of AZD8055, BEZ235, or GDC0941 using a HP D300 Digital Dispenser (Hewlett-Packard), while all other experimental conditions remained the same. Each assay was carried out in biological triplicate. Each replicate of a dose–response experiment was further analyzed by normalization to the negative and positive control (the normalized data are provided in Supplementary Table S3) and fitting to a four-parameter sigmoid function that allowed for the calculation of the IC₅₀ (dose at which viability is 50% of the untreated control). The IC₅₀ estimates are provided in Supplementary Table S4. For model inference, full dose–response curve data were used.

HCC1806 parental, GFP- and 4E-BP1-expressing cells were seeded at 600 cells/well, whereas the HCC1937 panel was seeded at 1,200 cells/well, in 96-well plates. Cells were treated, 24 hours after seeding, with a 9-point 1:3 dilution series of AZD8055 (top dose 3.3 × 10⁻⁶ mol/L) or BEZ235 (1.1 × 10⁻⁶ mol/L) using a HP D300 Digital Dispenser (Hewlett-Packard). Each condition, including an untreated negative control and a phenyl arsine oxide (1 × 10⁻⁶ mol/L) treated positive control, were set up in technical duplicate. Media and drugs were changed every 3 to 4 days over a period of 10 to 11 days of treatment. Cells were then washed with PBS, fixed with 3.7% formaldehyde/PBS, and stained in 0.1% crystal violet solution. Images of dried, stained cells were digitized on a Perfection V750 PRO scanner (Epson).

Steady-state RNA and protein expression was determined from cells seeded in 60-mm dishes and grown for 48 hours (seeding densities are listed in Supplementary Table S1). RNA expression was determined using RNA sequencing by the NKI Genomics Facility and protein expression was determined in biological triplicate using reverse-phase protein array (RPPA) analysis by the MD Anderson Cancer Center RPPA Facility. DNA was extracted from cell pellets and targeted sequencing was done using the Human Kinome Capture Kit (Agilent Technologies), containing baits for kinases and cancer-related genes. RNA-sequencing data are available at ArrayExpress, reference E-MTAB-4801, and the normalized read counts are provided in Supplementary Table S5. DNA sequencing data are available at the European Nucleotide Archive, reference PRJEB14120, and the mutation calls and estimated copy number profiles are provided in Supplementary Tables S6 and S7. The RPPA data are included as Supplementary Table S8.

Cells were seeded in 100-mm dishes (BT549 and CAL-120 at 2.5 × 10⁵; Hs 578T at 3 × 10⁵; HCC1806 at 4 × 10⁵; HCC1937 at 6.25 × 10⁵) and cultured for 48 hours, then treated with AZD8055 (1.11 × 10⁻⁷ mol/L), BEZ235 (3.7 × 10⁻⁸ mol/L), or vehicle (DMSO) for a further 24 hours. Cells were washed once with ice-cold PBS and lysed in lysis buffer (25 mmol/L Tris-HCl, pH 7.6, 1% Triton X-100, 1 mmol/L DTT) supplemented with cOmplete protease and phosSTOP phosphatase inhibitor cocktails (Roche). Lysates were cleared and assayed for protein concentration, then total protein samples were prepared using 20 μg of protein lysate. Cap pull-down samples were prepared by combining 50 μg of total lysate with 20 μL prewashed m⁷GTP-agarose (Jena Bioscience), made up to a total volume of 500 μL with lysis buffer and tumbled at 4°C overnight. The following day, cap pull downs were washed three times in ice-cold lysis buffer, then heated at 70°C for 10 minutes in 20 μL 1× Novex LDS Sample Buffer and sample reducing agent. The eluate from the cap pull downs and the total protein control samples were then immediately separated on Novex 4% to 12% gradient gels and immunoblotted using primary antibodies to 4E-BP1, eIF4G, and HSP90 (for total lysate samples only), then reprobed to detect eIF4E protein.

pLX304-4E-BP1 was obtained from the CCSB-Broad Lentiviral Expression Collection, whereas the pLX304-GFP control construct was generated as outlined previously (18). To produce lentiviral particles, HEK293T cells were cotransfected with the pLX304-4E-BP1 or -GFP bearing construct and a lentiviral packaging mix (pRSV-Rev, pMDLg/pRRE, pCMV-VSV-G; Addgene) using Polyethylenimine (PEI, Linear MW 25,000; Polysciences Inc.). Media were changed 24 hours after transfection. After a further 24 hours, viral supernatant was collected and 0.45-μm filtered. HCC1806 and HCC1937 cells were transduced in the presence of hexadimethrine bromide (Sigma-Aldrich) and after 48 hours selected using blasticidin.

HCC1806 parental, GFP-, and 4E-BP1–expressing cells were seeded at 800 cells/well, whereas the HCC1937 panel was seeded at 1,000 cells/well, in 384-well plates with 4 to 6 replicates per condition. Proliferation was monitored using the IncuCyte ZOOM instrument (Essen Biosciences).

An overview of the computational model and modeling procedure is given in the Quick Guide to Equations and Assumptions, and visualized in Fig. 1. The model inference was done using BCM (19). All files required for running the inference are included in Supplementary File S1. A detailed description of all equations, data preprocessing, and inference algorithms is given in the Supplementary Materials and Methods.

We set out to establish an integrative computational model capable of explaining observed therapeutic responses based on molecular measurements. To this end, we sourced and comprehensively characterized 30 breast cancer cell lines (Fig. 2; Supplementary Table S1). Given the need for targeted treatment options for the triple-negative breast cancer subtype, the panel was enriched for triple-negative cell lines (18), with four ER⁺, four HER2⁺, and four ER⁺/HER2⁺ cell lines included to represent the other major subtypes.

The panel was characterized for response to seven kinase inhibitors, including AZD8055 (mTOR inhibitor), BEZ235 (dual mTOR/PI3K inhibitor), GDC0941 (PI3K inhibitor), MK2206 (AKT inhibitor), PD0325901 (MEK inhibitor), lapatinib (dual EGFR/HER2 inhibitor), and foretinib (cMET/VEGFR2 inhibitor). The sensitivity of each cell line to these inhibitors was determined in 10-point, 72-hour dose–response assays, in biological triplicate and technical quadruplicate (summarized by IC₅₀ values in Fig. 2; Supplementary Fig. S2). The drug sensitivity measurements largely agree with those obtained in the Genomics of Drug Sensitivity in Cancer screen (Supplementary Fig. S3; refs. 10, 17, 20). However, our use of a more focused panel of 30 cell lines and seven drugs allowed us to obtain more precise measurements (Supplementary Fig. S4).

In addition to response data, we profiled the panel for mutation and copy number by DNA-seq, RNA expression by RNA-seq, protein expression, and phosphorylation by RPPA as well as proliferation rate under untreated, steady-state growth conditions (Fig. 2, bottom right; see Supplementary Figs. S5–S9, for enlarged versions of the graphs). This molecular characterization was done in the absence of drug treatment. The cell lines harbor a range of genetic events that occur in breast tumors and are present at comparable frequencies (see Supplementary Table S13, for a comparison of mutation frequencies with tumors from The Cancer Genome Atlas; ref. 2). This cell line panel thus represents a relevant model system for the genetic diversity in breast tumors.

To first understand which signaling is relevant for each drug, we developed a modeling framework, inference of signaling activity (ISA), to infer the signal strengths and signaling activities from all available data, as described in the Quick Guide to Equations and Assumptions and further detailed in the Supplementary Materials and Methods. We constructed a literature-based model and first fitted this to the response data for each drug separately, in conjunction with the molecular data measured in untreated cells. Although all of the interactions included in the model are well documented (Supplementary Tables S9 and S10), their relative contribution or significance is not known. For example, activating mutations in PIK3CA, the loss of PTEN, or the expression of growth factors can all lead to activation of the PI3K pathway. However, it is unclear whether their effects are equally important, and if not, which of them has a stronger effect in a particular context.

Figure 3A illustrates the model estimates of signaling strengths (the links between signaling molecules) for lapatinib treatment. Values of the strength parameters indicate which signaling connections are important for propagating an oncogenic signal down to the proliferation node. For example, the link between ERBB2 amplification and ERBB2 activation has a strong peak at nonzero values (the density plot of ERBB2amp→ERBB2), thus indicating that the ERBB2 amplification gives rise to a proliferative signal. It is well known that amplification of ERBB2 and the resulting overexpression and auto-activation of this receptor provides a strong proliferation signal (21), and that this signal can be inhibited by lapatinib (22). The model provides estimates for the downstream signaling (e.g., indicating that ERBB2 signals more to PI3K than to CRAF) and for the contribution of each of the resistance mechanisms. PIK3CA mutations indeed contribute to the proliferative signal (PIK3CAhelical→PI3K), and the model predicts that PIK3R1 mutations may have a similar effect, as we see that the parameter describing how much a PIK3R1 mutation activates PI3K tends toward higher values, and is most likely nonzero (PIK3R1→PI3K). However, the uncertainty in this parameter is large, because there is only one cell line that carries such a mutation, and consequently this parameter is only weakly constrained. As a last example, the contribution of HGF autocrine signaling (9) is represented by the parameter controlling how strongly expression of HGF leads to activation of the MET receptor (HGF→MET). The posterior probability distribution of this parameter closely follows the prior, indicating that this parameter, and thus the importance of this potential resistance mechanism, cannot be determined from the current data. Together, this shows that the ISA modeling approach can be used to infer the contribution of different components driving sensitivity and resistance from all available data, while taking into account whether the parameters are identifiable.

We can use the estimates of the signaling activities (values of the nodes) to further explore the difference in signaling flow and drug response between cell lines. Figure 3B shows the estimates for ERBB2 activity and PIP3 activity in the lapatinib-treated condition scattered against the untreated condition. From the left panel, it can be seen that only the eight ERBB2-amplified cell lines show ERBB2 signaling activity, and that this activity is reduced upon lapatinib treatment in all these cell lines. In the right panel, we can see that in the lapatinib-sensitive cell lines, especially the most sensitive ones BT-474, SK-BR-3, and ZR-75-30, the reduction in ERBB2 activity also leads to a strong reduction in the PIP3 signal, whereas in the other cell lines, PIP3 signal persists, especially for the lapatinib-resistant cell line HCC1569 (see Supplementary Fig. S2, for the drug sensitivity estimates). Two non-ERBB2-amplified cell lines also have a reduced PIP3 signal upon lapatinib treatment, including T-47-D and HCC1806, which stems from their inferred EGFR activity. This illustrates the utility of the model, given that there were no molecular measurements collected in the treated conditions, and these signaling estimates are inferred from the untreated molecular data combined with the relative viability data in the treated condition. A comparison of the inferred model estimates with molecular measurements in treated condition is described later, after all model adaptations have been considered.

As a second example, for the mTOR inhibitor AZD8055, we find several factors that are associated with response. As expected, PIK3CA mutations are strongly activating in this context and cell lines are apparently dependent on this activation (Fig. 3C, PIK3CAkinase/helical→PI3K and mTORC1→proliferation), which has previously been shown (23). In addition, we find that MYC activation, as a result of gene amplification, can provide a resistance mechanism to this mTOR inhibitor (MYCamp→proliferation), providing another validation of our approach (24).

To facilitate the further exploration of all the signaling estimates, we generated an interface that is available at http://ccb.nki.nl/software/BCCL_KI_response_model/. This website displays the signaling strengths in each cell line upon exposure to drugs. Figure 3D shows an example of BT-474 treated with lapatinib. In this case, the model indicates that, for example, the MAPK pathway is barely involved, that there is no drug resistance provided by this cell line's RPS6KB1 amplification, and that instead, lapatinib mainly inhibits the PI3K pathway.

Although the signaling estimates appear to be reasonable, it is useful to have a systematic test of how well the fitted model describes the data. For this we used the posterior predictive distribution, which describes a new, predicted dataset based on the fitted model. We can overlay this predicted dataset on the observed measurements to have a convenient way of identifying which measurements can and cannot be explained by the model. Note that the posterior predictive distribution is not used as a measure of out-of-sample prediction; that is, it does not test how well the model predicts the behavior of unseen cell lines. Instead, it is used as a measure of goodness of fit, allowing an exploration of whether the model can describe the behavior of the cell lines at hand.

Fig. 4A shows the posterior predictive distribution overlaid on the measurement data for lapatinib. It is clear that the present model can accurately describe the relative proliferation of the cell lines as a function of drug concentration for almost all cell lines. For example, the sensitivity of the ERBB2-amplified line BT-474 and the resistance of the PIK3CA-mutated and ERBB2-amplified line MDA-MB-361 (MM361) can both be recapitulated by the model. Overviews of all posterior predictive checking for the drug response and phosphorylation data are supplied in Supplementary Data S1.

Although the model explained most of the drug response variability for lapatinib, we noticed that for some drug–cell line combinations, the fit was not as precise. For example, for foretinib, an inhibitor of c-Met and VEGFR2, we noticed that one cell line in particular, MFM-223, was much more sensitive than the model could describe (Fig. 4B). We therefore investigated the experimental data to find a possible reason for this discrepancy. A discrepancy in a single cell line is not sufficient to apply statistical tests, but we did note that this cell line has a strong FGFR2 amplification. We therefore searched the literature to see whether there is a connection and found that foretinib has in fact been reported to inhibit FGFR2 in addition to its original design targets (25). When this additional target of foretinib is added to the model, we indeed obtain a significantly improved fit (Fig. 4B). The sensitivity of other cell lines like CAL-51 and HCC-1187 to foretinib is still not explained exactly, but we have not found other potential explanations for this in the data or literature, and therefore further studies would be needed to address this.

We also noticed that the model was not able to explain the response of some cell lines to mTOR inhibitors (Fig. 4C), but were unable to find additional mechanisms in the literature that could explain these discrepancies. Several cell lines are sensitive to these inhibitors even though they do not possess any of the factors known to cause sensitivity, and conversely some cell lines are resistant despite having such sensitizing factors. For example, although BT-549, HCC1395, and HCC1937 have all lost PTEN expression, only BT-549 is sensitive to BEZ235 treatment.

We therefore further interrogated the dataset to find additional drug sensitivity mechanisms with which we could extend the model to better explain the experimental observations. With multiple sensitive cell lines we can use statistical tests. We divided the cell line panel in groups of sensitive and resistant lines using Gaussian mixture modeling, and tested whether any genes were differentially expressed between these groups at either the RNA or protein level. The full lists of differentially expressed genes for all drugs are given in Supplementary Table S14. To further filter the differentially expressed genes, we also calculated their distance to the signaling molecules included in the model using protein–protein interaction networks (26). This provided potential candidate regulators that are not only differentially expressed, but are also functionally closely related to the signaling molecules in the model (listed in Supplementary Table S14). For both mTOR inhibitors (AZD8055 and BEZ235), the protein expression level of 4E-BP1, in addition to 4E-BP1 phosphorylation level, showed the strongest differential expression (Fig. 4D). At the RNA level, differential expression of EIF4EBP1 was also associated with BEZ235 response. Together, these data indicated that cell lines with high expression of this protein are more sensitive to mTOR inhibitors (Supplementary Fig. S10).

To test whether the inclusion of this factor provides a better explanation for the sensitivity of some of the lines to mTOR inhibitors, we expanded the model to include the protein expression levels of 4E-BP1. Although 4E-BP1 as a node was already included as a downstream target of mTORC1 in our previous models, only its phosphorylation state, not its protein expression level, was taken into account. Specifically, in Eq. A, the variable E_i,j was previously based only on the binarized RNAseq expression data (Supplementary Fig. S11), and we modified this to include the RPPA protein expression levels (Supplementary Materials and Methods for details). Fig. 4E shows the model with the protein expression levels of 4E-BP1 included, whereas Fig. 4C shows the results without 4E-BP1 included. The posterior predictive checking (bottom panel) clearly shows that the expanded model provides an improved fit to the data for multiple cell lines (especially for CAL-120, BT-549, CAMA-1, and ZR-75-30), while not compromising the fit of other cell lines. The log Bayes factor between the two models (the difference in log evidence) is 226 in favor of the expanded model, where a log Bayes factor greater than 5 indicates very strong evidence (27). This indicates that the difference is highly significant and that the improvement in fit is not merely the result of adding more free parameters. Fig. 4F shows the posterior probability of two parameters of 4E-BP1 signaling in the expanded model. The 4E-BP1 expression coefficient, shown in the top panel, describes how strongly the protein expression level of 4E-BP1 affects the amount of signal transmitted. Because the value of this parameter is very high, it indicates that the total protein expression level is a strong limiting factor for 4E-BP1 signaling in response to mTOR inhibitors. The bottom panel shows the parameter controlling the strength of 4E-BP1 signaling to proliferation, and as this is also found to be nonzero with high certainty, it implies that the 4E-BP1 signal is important for determining proliferation rate under mTOR inhibitor treatment. The computational analysis therefore predicts that the protein expression level of 4E-BP1 is an important factor in explaining mTOR inhibitor sensitivity, in addition to the already known factors determining sensitivity and resistance.

The model uses pretreatment molecular data and measurements of relative viability after drug treatment to infer signaling activities after drug treatment. To gain more confidence in the signaling estimates, we compared the estimates with measurements of protein phosphorylation of cells while on treatment (28), and found that they generally agree (Supplementary Note S1; Supplementary Fig. S12). We also constructed a reduced model, which could be fitted to the response data of all seven inhibitors at the same time (Supplementary Note S2; Supplementary Fig. S13). Finally, using an extended version of the modeling formalism (29), we confirmed that even though feedback signaling is likely to be active, the inclusion of such feedback loops does not affect how well the variability in drug response can be described (Supplementary Note S3; Supplementary Fig. S14).

Intrigued by the model prediction that 4E-BP1 protein expression is associated with mTOR inhibitor response, we investigated the biological effect of 4E-BP1 expression directly. For this, we turned to a subset of our panel of breast cancer cell lines that showed differences in response to AZD8055 and BEZ235 (Fig. 5A). These included three of the most mTOR inhibitor sensitive cell lines (BT-549, CAL-120 and Hs 578T) all bearing a gain in the EIF4EBP1 gene-containing genomic region (Supplementary Fig. S15), which also express high levels of 4E-BP1 protein (Fig. 5B; Supplementary Fig. S16), and two insensitive cell lines (HCC1806 and HCC1937) that do not harbor a gain of the EIF4EBP1 locus and express low levels of 4E-BP1 protein. Given that high 4E-BP1 expression may drive cells to recalibrate signaling in the pathway by increasing the expression and/or activity of mTOR, we investigated this possibility further. We first checked whether expression of 4E-BP1 and mTOR were correlated (Supplementary Fig. S17A). Although three of the most highly 4E-BP1–expressing cell lines do show an increase in mTOR expression at the protein level, in the lines chosen for our functional studies, only CAL-120 shows elevated mTOR expression. The remaining four lines (BT549, HS578T, HCC1806, HCC1937) show comparable expression of mTOR, despite vast differences in both 4E-BP1 expression and response to mTOR inhibitors. At the RNA level (Supplementary Fig. S17B), these associations were lost, suggesting that a concurrent posttranscriptional upregulation of mTOR is not a general mechanism of mTOR inhibitor sensitivity in 4E-BP1–overexpressing cells. We then investigated mTOR signaling in more detail by analyzing the phosphorylation and protein expression levels of several members of the PI3K, mTOR, and MAPK pathways following 24-hour treatment with two mTOR inhibitors, AZD8055 and BEZ235, in these five cell lines. This showed that both compounds had effective on-target activity, leading to reduced phosphorylation of AKT (S473), S6 (S235/236), and 4E-BP1 (S65) across all five lines, at similar compound concentrations (Fig. 5B). A minor compensatory increase in ERK phosphorylation was detected following inhibitor treatment. These data suggest that the difference in mTOR inhibitor sensitivity between these five cell lines is not caused by a difference in the compounds' ability to inhibit mTOR signaling in these lines.

To uncover the mechanism via which 4E-BP1 protein expression levels could affect response to mTOR inhibitors, we investigated the effect of inhibitor treatment on the formation of the eIF4F translation initiation complex (extensively reviewed in refs. 30, 31). As illustrated in the schematic in Supplementary Fig. S18A, the eIF4F translation complex is composed of the eIF4E and eIF4G proteins, among others. 4E-BP1 is known to negatively regulate this complex by binding and sequestering the eIF4E subunit. This displaces eIF4G from binding to eIF4E and as a result the eIF4F complex cannot initiate cap-dependent translation. The sequestering of eIF4E by 4E-BP1 is, however, inhibited when 4E-BP1 is phosphorylated by mTORC1 on several sites, which is the case when nutrients and growth factors are not limiting. Under nutrient or growth factor depletion, or alternatively following treatment with mTOR inhibitors, 4E-BP1 becomes dephosphorylated, binds to eIF4E, and thus eIF4F complex activity is repressed (Supplementary Fig. S18B). This leads to the inhibition of translation, most acutely for mRNAs with complex 5′ untranslated regions (UTR) that include proliferation, survival, and tumor-promoting genes.

We investigated the dynamics of these interactions in the three mTOR inhibitor sensitive and two insensitive cell lines using the m⁷G-cap pull-down assay, which allows the visualization of the changes in eIF4G or 4E-BP1 binding to eIF4E following treatment, as compared with their expression in the total protein lysate (Fig. 5C). Our results show that mTOR inhibitor treatment leads to an increase in the binding of 4E-BP1 to eIF4E in all five cell lines, irrespective of their mTOR inhibitor response profile. In the sensitive cell lines, this increase in 4E-BP1 binding was sufficient to decrease eIF4G binding to eIF4E, as expected. In the insensitive cell lines though, the binding of 4E-BP1 to eIF4E was unable to displace eIF4G from eIF4E. This suggests that 4E-BP1 protein expression in the insensitive cell lines is below a critical threshold needed to effectively inhibit eIF4F complex formation following mTOR inhibitor treatment (Supplementary Fig. S18C), and likely explains the difference in mTOR sensitivity between these two sets of cell lines.

To further investigate whether an increase in 4E-BP1 protein expression is sufficient to increase mTOR inhibitor response, we used a lentiviral vector to overexpress 4E-BP1 in the two insensitive cell lines, HCC1806 and HCC1937 (Fig. 6A). The 4E-BP1 protein, as well as a GFP control, was effectively overexpressed in both lines and the former was detectably phosphorylated. The expression of either protein did not affect the proliferation of these cell lines, suggesting that the activity of overexpressed 4E-BP1 was efficiently inhibited by mTORC1-mediated phosphorylation (Fig. 6B). We next tested the impact of increased 4E-BP1 protein expression on the sensitivity of the cell lines to the mTOR inhibitor AZD8055, the dual mTOR/PI3K inhibitor BEZ235, as well as the PI3K inhibitor GDC0941. As shown in Fig. 6C, in short-term 72-hour drug treatment assays, the 4E-BP1-overexpressing cell lines were markedly more sensitive to mTOR inhibitors, responding at lower drug concentrations and with a decreased overall survival at higher drug concentrations. In contrast to the mTOR inhibitors, sensitivity of the 4E-BP1-overexpressing cell lines to the PI3K inhibitor GDC0941 was not increased, implying that it is specifically the inhibition of mTOR activity that is beneficial in improving response of highly expressing 4E-BP1 cell lines. We were also able to validate that 4E-BP1 overexpression increased sensitivity to mTOR inhibitors over a longer treatment period, namely 10 days, as shown in Fig. 6D.

Together, the above results show that the level of 4E-BP1 protein expression is a determinant of sensitivity to mTOR inhibitors in breast cancer cell lines, and illustrates the utility of our computational model in identifying novel determinants of drug response in cell lines.

Cell line panels have the potential to provide us with a better understanding of the variability in drug response between patients. Previous efforts of linking molecular characteristics to drug sensitivity in cell line panels have identified several known and novel associations (10, 11, 13–16). Here we showed that by combining extensive measurements with mathematical modeling, a more detailed understanding of variability in drug sensitivity can be achieved. The use of Bayesian statistical analysis allowed for the simultaneous integration of diverse data types with prior knowledge from the literature.

Models are by definition a simplified representation of the system. Two major simplifications that were used here are the assumption of quasi-steady state and the absence of feedback signaling. The benefit of using these simplifications is that significantly more components of cellular signaling can be included in the model. It is reassuring that a model with these simplifications can describe a large part of the variability in short-term drug response. Studying longer-term drug response and, for example, adaptive resistance will likely require the incorporation of dynamics and feedback mechanisms. Using dynamical models or the inclusion of feedback mechanisms does not pose any theoretical problems for the modeling approach we used. However, the computational cost is significantly higher, and additional intervention or time-course data would be needed to constrain the parameters.

Recently, Fey and colleagues described a computational model of JNK signaling, containing five signaling proteins that provided prognostic information in patients with neuroblastoma (32), showing that such computational models may be useful also in a clinical setting. Their model was informative in a specific subset of patients, namely those with MYCN-amplified tumors, whereas in the general population individual biomarkers were still more informative. This indicates that it is necessary for models to incorporate various different oncogenic drivers, residing in multiple signaling pathways, to capture the variability across a wide range of patients.

A method that integrates multiple data types to obtain pathway activation status is PARADIGM (33). Heiser and colleagues have used a modified version of this method, SuperPathway analysis, to link pathway activation status to drug response in a large breast cancer cell line panel (12). They found, among others, that upregulation of DNA damage response pathways was associated with sensitivity to cisplatin, although this was not further tested experimentally. Indeed, PARADIGM and SuperPathway analyses do not shed light on how this association might work mechanistically, and the involvement of individual components of the DNA damage response pathway, such as TP53, ATM, and BRCA1/2, is not investigated. In contrast, the approach presented here provides detailed signaling flows and estimates the relative contributions from all drivers and sensitivity mechanisms included in the model, making the in silico findings amenable to experimental validation.

In the course of refining our model, we found that elevated 4E-BP1 protein expression played an important role in the response of breast cancer cell lines to AZD8055 and BEZ235, two inhibitors targeting the mTOR kinase. We further validated this observation in 4E-BP1 overexpression studies, showing that this provides a pool of an endogenous translation inhibitor available for activation, and thus inhibition of cap-dependent translation, via mTOR inhibitor treatment.

On a mechanistic level, these findings are in line with prior work using transformed mouse embryonic fibroblasts, which showed that the ratio of eIF4E/4E-BP1 expression can predict response to mTOR-directed therapy (34). That is, a higher ratio of eIF4E/4E-BP1 predicted poorer response to mTOR inhibitors. Consistent with this, EIF4E amplification was reported as a mechanism of AZD8055 resistance in a SW620 colorectal cell line model (35). Conversely, a lack of 4E-BP1 expression in lymphoma cells, thus a reduced ability to restrain eIF4E activity, has been shown to lead to resistance to mTOR inhibition, an effect reversed by exogenous expression of the 4E-BP1 protein (36), similar to our findings in 4E-BP1-overexpressing HCC1806 and HCC1937 cell lines. An exception to these findings is the report of elevated 4E-BP1 protein levels in mTOR inhibitor treatment-resistant luminal subpopulation of prostate cancer cells (37), although it remains to be investigated whether this observation is restricted to prostate cancer or the luminal subtype. Most recently, a study by Wang and colleagues (38) added further weight to the hypothesis that elevated 4E-BP1 expression can be a marker of mTOR inhibitor sensitivity. They show that the combination of an mTOR inhibitor and an HDAC inhibitor, the latter acting to derepress Snail-mediated 4E-BP1 transcriptional inhibition, can synergize to inhibit tumor growth in mice.

Notably, in our cell line panel, the overexpression of 4E-BP1 resulted from a copy number gain in the genomic region encoding EIF4EBP. Focal amplification of the 8p11-12 region, the region containing EIF4EBP1, is a known event in breast cancer occurring in approximately 15% of cases, and patients who harbor this event in their primary tumor have a much higher likelihood of relapse (39). Previous studies have identified various genes in this region as potential oncogenes, with most evidence so far supporting FGFR1 (40–42) and ZNF703 (43–45). These two genes lie close to and on either side of EIF4EBP1. Our cell line experiments show that overexpression of 4E-BP1 alone does not affect viability in vitro. Although we cannot exclude that amplification of 4E-BP1 can contribute to a transformed phenotype in tumors (such as affecting invasion, migration, or cell viability in vivo), various studies indicate that the cellular function of 4E-BP1 is consistent with a role as a tumor suppressor gene, rather than as an oncogene (34–36, 46–49). It therefore seems plausible that the amplification of EIF4EBP1 is a passenger event. This raises the possibility that by selecting for amplification of nearby oncogenes, the tumors have also introduced a specific “passenger vulnerability”—that is, a passenger aberration that introduces a vulnerability to a particular drug. If the drug sensitivity association translates from cell lines to patients, testing for EIF4EBP1 amplifications could identify patients who may benefit most from treatment with mTOR inhibitors such as AZD8055 and BEZ235.

One difference between the cell line and patient data is that the subtypes in which the gain of EIF4EBP1 is present vary. Although we identified this event predominantly in the triple-negative and ER+ breast cancer cell lines, in patients, the 8p11-12 amplicon is present almost exclusively in ER⁺ tumors. Interestingly, phase III clinical trials with the allosteric mTOR inhibitor, everolimus, have been carried out in the ER⁺ setting, showing a significant improvement in progression-free survival (PFS) in patients who received everolimus versus placebo in addition to the aromatase inhibitor, exemestane (50). A subsequent analysis of this data, exploring associations of PFS with common genetic aberrations, has found no improvement in PFS in patients with FGFR1 gene amplification (generally coamplifying the EIF4EBP1 gene; ref. 51). Although this may suggest that there is no increased clinical benefit for mTOR inhibitors in patients with 8p11-12 amplified tumors, the inhibitory activity of everolimus, a rapalog, versus active-site inhibitors, such as AZD8055 and BEZ235, differs substantially. Although everolimus is an allosteric inhibitor of predominantly mTOR complex 1 (mTORC1), active-site inhibitors are able to target both mTORC1 and mTORC2 (30). Most importantly, the active-site inhibitors have been shown to result in a much more potent inhibition of downstream mTOR signaling, specifically with respect to the inhibition of 4E-BP1 phosphorylation, whereas in the case of rapalogs this phosphorylation is restored within hours of treatment (52, 53). As such, it is likely that the absence of an association between amplification in the 8p11-12 region containing the FGRF1 and EIF4EBP1 genes and PFS with everolimus treatment observed by Hortobagyi and colleagues (51) stems from limited inhibition of 4E-BP1 phosphorylation. Together, these studies suggest that the benefit of rapamycin treatment in ER⁺ tumors results from a mechanism distinct to that of 4E-BP1 inhibition, but also emphasize the need for further study to determine the efficacy of active-site mTOR inhibitors in patients with 8p11-12 amplifications.

We conclude that the combination of mathematical modeling of signaling in response to drug treatment with a large panel of molecularly characterized cell lines demonstrates the usefulness of combining large data sets with prior knowledge to uncover key determinants of drug sensitivity. Although the work presented here applied the ISA methodology to explain drug response in a panel of cell lines, our ultimate goal is to develop this approach into a tool for predicting response in patients. Further validation on clinical samples is of course required, but a clear benefit is that the molecular characterization, needed as input for a model, can be performed directly on biopsy material. This circumvents a need for ex vivo culture or lengthy response profiling, making the approach amenable to clinical application. We believe that such systematic, quantitative approaches to the understanding of drug responses hold promise as tools for achieving the goal of precision medicine in cancer.

No potential conflicts of interest were disclosed.

Conception and design: K. Jastrzebski, B. Thijssen, R.L. Beijersbergen, L.F.A. Wessels

Development of methodology: K. Jastrzebski, B. Thijssen

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): K. Jastrzebski, K. de Lint, I.J. Majewski

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): K. Jastrzebski, B. Thijssen, R.J.C. Kluin, L.F.A. Wessels

Writing, review, and/or revision of the manuscript: K. Jastrzebski, B. Thijssen, K. de Lint, R.L. Beijersbergen, L.F.A. Wessels

Study supervision: R.L. Beijersbergen, L.F.A. Wessels

We are grateful to Thomas Kuilman for performing the CopywriteR analysis and the NKI Genomics Facility for performing high-throughput sequencing experiments and subsequent data processing. We are grateful to Jordi Vidal Rodriguez for help with the drug response assays and maintenance of the cell line panel. We thank the Beijersbergen, Wessels, and Bernards groups for helpful discussions. This work was funded by the Cancer Systems Biology Center grant (853.00.120) from the Netherlands Organization for Scientific Research (NWO) and the Cancer Genomics Netherlands Program supported by the Gravitation program of the Netherlands Organization for Scientific Research (NWO).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.