Mathematical modelling of breast cancer cells in response to endocrine therapy and Cdk4/6 inhibition.

Oestrogen receptor (ER)-positive breast cancer is responsive to a number of targeted therapies used clinically. Unfortunately, the continuous application of any targeted therapy often results in resistance to the therapy. Our ultimate goal is to use mathematical modelling to optimize alternating therapies that not only decrease proliferation but also stave off resistance. Toward this end, we measured levels of key proteins and proliferation over a 7-day time course in ER+ MCF-7 breast cancer cells. Treatments included endocrine therapy, either oestrogen deprivation, which mimics the effects of an aromatase inhibitor, or fulvestrant, an ER degrader. These data were used to calibrate a mathematical model based on key interactions between ER signalling and the cell cycle. We show that the calibrated model is capable of predicting the combination treatment of fulvestrant and oestrogen deprivation. Further, we show that we can add a new drug, palbociclib, to the model by measuring only two key proteins, cMyc and hyperphosphorylated RB1, and adjusting only parameters associated with the drug. The model is then able to predict the combination treatment of oestrogen deprivation and palbociclib. We illustrate the model's potential to explore protocols that limit proliferation and hold off resistance by not depending on any one therapy.

Introduction

Despite advances in treatment options, it is estimated that 42 260 women and men will die as a result of breast cancer in the USA this year [1]. Seventy per cent of all breast cancer cases are oestrogen receptor positive (ER+) and targeted therapies such as endocrine therapy or Cdk4/6 inhibitors are used to treat this clinical subtype, resulting in dramatic improvements in long-term survival rates. Endocrine therapies act to decrease ER signalling in a variety of ways: (i) by depriving the ER of its ligand oestrogen (17β-oestradiol, hereafter referred to as E2) via aromatase inhibitors (e.g. letrozole), (ii) by competitively inhibiting the binding of E2 to the ER (the mode of action of tamoxifen), or (iii) by increasing the degradation of the ER via the proteasome (the mode of action of ICI 182 780 (ICI; Faslodex/fulvestrant)). This inhibition of ER signalling can halt proliferation by arresting cells in G1 and inducing cell death [2-4]. Inhibitors of cyclin-dependent kinases Cdk4/6, such as palbociclib, have been used in conjunction with endocrine therapies to provide a more durable inhibition of proliferation [5-7].

Unfortunately, long-term application of a targeted therapy may lead to the development of resistance and recurrence of cancer [8-12]. However, sequential application of various therapies and drug holidays, in a repeating cycle, may provide a means of delaying resistance. To develop such a sequential therapy, in light of the many possibilities for drugs, doses and timings, mathematical modelling of the system is likely to be critical [13]. Ultimately, such a model should be able to incorporate numerous drugs, predict the proliferation of ER+ cells in response to these drugs over relatively long time scales and account for some of the many possible resistance mechanisms. Such a model would allow for optimizing sequential therapies, subject to numerous constraints, to propose therapeutic protocols that can be evaluated experimentally.

Our study takes the first step along this road by building a model, based on known mechanisms from the literature, to describe the changes in protein expression and proliferation of an E2-dependent, asynchronous population of ER+ MCF-7 breast cancer cells in response to various clinically relevant therapies over a time period of one week. The therapies we consider are E2 deprivation (a surrogate for aromatase inhibitors), ICI (a selective ER degrader) and palbociclib (a Cdk4/6 inhibitor). Because the model has a mechanistic basis, we show that it is relatively straightforward to add additional drugs relevant to the mechanism and to predict the effects of combinations of drugs. We also show how the model might be used to develop protocols that delay resistance, once it has been validated at longer time scales and the dynamics of resistance acquisition are better understood.

While there have been a number of previous mathematical models dealing with ER signalling and its impact on the cell cycle [14-16], these models differ from ours in that they are phenomenological and not experimentally calibrated or validated. There are also a number of mathematical models that describe the G1–S transition in single cells [17-21]. Most of these investigate the bistable switch governing the commitment to proliferation at the restriction point [17]. While our model is based on a similar structure, our aim is to capture the behaviour of an asynchronous population of cells, creating an average model, so we do not expect our parameterization to create a bistable switch or be valid for single cells.

Before presenting our results, it is important to clarify what can be expected of a model of the type we propose. The behaviour of protein levels in cells in response to therapeutic perturbations is the result of thousands, if not millions, of interactions among DNA, RNA and proteins. It is clearly impossible at present to model this level of complexity, and so of necessity we use a vastly simplified model to attempt to capture the key effects of therapy. This has two important implications for the modelling. First, many reaction rates used in the model account for numerous unmodelled effects (e.g. the rate at which ligand-bound ER upregulates the transcription factor c-Myc), and thus there is no corresponding rate in the literature. Second, the simple model is unlikely to provide an excellent match to the experimental data. Therefore, our goal is to produce a model that fits the data as well as possible and then show that the model is useful for its intended purpose: predicting the results of combination therapies and making it easy to add new drugs that operate on the mechanisms being modelled.

Material and methods

Cell culture and reagents

MCF-7 cells were obtained from Tissue Culture Shared Resources at Lombardi Comprehensive Cancer Center, Georgetown University, Washington, DC [22,23]. MCF-7 cells depend on E2 for growth, and therefore, to enable us to control the level of E2 in the cell culture medium, cells were grown in phenol red-free improved minimal essential medium (Life Technologies, Grand Island, NY; A10488-01) with 10% charcoal-stripped calf serum (CCS) and supplemented with 10 nM 17β-oestradiol (oestrogen; E2). E2 deprivation was obtained by washing cells 24 h post-plating (t = 0) with phosphate-buffered saline (PBS) and adding complete medium without E2 for the indicated times. ICI (Faslodex/Fulvestrant; ICI182,780) and palbociclib were obtained from Tocris Bioscience (Ellisville, MO). MCF-7 cells were authenticated by DNA fingerprinting and tested regularly for Mycoplasma infection. All other chemicals were purchased from Sigma-Aldrich (St. Louis, MO).

Cell proliferation assays

Cells were seeded at a density of 1–2 × 105 cells well−1 or 4–5 × 104 cells well−1 in 100 mm or 60 mm plates, respectively. Cells were then trypsinized, resuspended in PBS and counted using a Z1 Single Coulter Counter (Beckman Coulter, Miami, FL). Three independent experiments were done. Data were normalized to cell number at t = 0 and are presented as the mean ± s.e. from all three experiments.

Western blot analysis

For Western blot analysis, cells were lysed for 30 min at 4°C in lysis buffer (50 mM Tris-HCl, pH 7.5, containing 150 mM NaCl, 1 mM EDTA, 0.5% sodium deoxycholate, 1% IGEPAL CA-630, 0.1% sodium dodecyl sulfate (SDS), 1 mM Na3VO4, 44 µg ml−1 phenylmethylsulfonyl fluoride) supplemented with Complete Mini protease inhibitor mixture tablets (Roche Applied Science). Total protein was quantified using the bicinchoninic acid assay (Pierce). Whole-cell lysate (20 µg) was resolved by SDS–polyacrylaminde gel electrophoresis. The following antibodies were used: c-Myc (no. 5605), cyclinD1 (no. 2978), p21 (no. 2947) and RB1 (no. 9309) were from Cell Signaling, Danvers, MA; ESRα (MA5-14104) was from Invitrogen; RB1-phosphorylated on Ser612 (OAAB16108) was from Aviva Systems Biology, San Diego, CA; actin (sc-47778) was from Santa Cruz Biotechnology, Santa Cruz, CA; β-tubulin (T7816) was from Sigma, St. Louis, MO.

Data collection

Table 1 summarizes the data that were collected. Protein levels were analysed by western blotting and cell numbers were determined using a Coulter counter. Protein levels and cell numbers were first normalized to 0 h. Protein levels in treatment conditions were then normalized to the +E2 case (control; untreated, grown with E2 in medium). See the electronic supplementary material for details.

Table 1.

Time points at which the various protein species and cell numbers were measured in response to the various treatments. +E2 (untreated, grown with E2 in medium), +E2 + ICI (treated with ICI in +E2 medium), –E2 (deprived of E2), –E2 + ICI (treated with ICI in –E2 medium), +E2 + palbociclib (treated with palbociclib in +E2 medium) and –E2 + palbociclib (treated with palbociclib in –E2 medium). Underlined time points were used for normalization; italicized for model parameter calibration; and bolded for model validation. h: hours and d: days.

condition	protein name	protein time points	cell number time points
+E2 (control)	ER, c-Myc, cyclinD1, p21, RB1-pp, total RB1	0 h, 4 h, 1d, 3d, 6d, 7d	0 h, 1d, 3d, 6d, 7d
+E2 + ICI	ER, c-Myc, cyclinD1, p21, RB1-pp, total RB1	0 h, 4 h, 1d, 3d, 6d, 7d	0 h, 1d, 3d, 6d, 7d
–E2	ER, c-Myc, cyclinD1, p21, RB1-pp, total RB1	0 h, 4 h, 1d, 3d, 6d, 7d	0 h, 1d, 3d, 6d, 7d
–E2 + ICI	ER, c-Myc, cyclinD1, p21, RB1-pp, total RB1	0 h, 4 h, 1d, 3d, 6d, 7d	0 h, 1d, 3d, 6d, 7d
+E2 + palbo	c-Myc, RB1-pp	0 h, 1d, 3d, 6d, 7d	0 h, 1d, 3d, 6d, 7d
–E2 + palbo	c-Myc, RB1-pp	0 h, 1d, 3d, 6d, 7d	0 h, 1d, 3d, 6d, 7d

Model structure

The interactions of key elements governing the biological mechanism associated with the effect of E2 signalling on proliferation are shown in figure 1a. E2 bound to ER enhances the production of cyclinD1 and c-Myc. CyclinD1 with its kinase partner Cdk4/6 drives RB1 to a hypophosphorylated state, RB1-p. CyclinE with its kinase partner Cdk2 drives RB1-p to a hyperphosphorylated state, RB1-pp, in which it can no longer bind and inhibit E2F, allowing E2F to drive the G1–S transition and proliferation. p21, representing both itself and p27, inactivates both cyclinD1:Cdk4/6 and cyclinE:Cdk2, but the increased production of c-Myc increases its inhibition of p21 production, helping to take the brakes off proliferation. Further details and references are provided in the electronic supplementary material.

Figure 1.

Wiring diagrams of the biological mechanism and the mathematical model. Solid lines with balls represent reversible binding reactions, other solid lines represent production or degradation of proteins, dashed lines represent influences with arrowheads representing enhancement and blunt heads representing inhibition. (a) Biological mechanism consisting of the following processes: (1) ER binds its ligand E2; (2) ICI binds ER and enhances its degradation; (3) ER:E2 enhances production of c-Myc; (4) ER:E2 enhances production of cyclinD1; (5) c-Myc inhibits the production of p21; (6) cyclinE binds to Cdk2; (7) p21 binds to cyclinE:Cdk2 and inactivates it; (8) cyclinD1 binds to Cdk4/6; (9) p21 binds to cyclinD1:Cdk4/6 and inactivates it; (10) palbociclib binds to Cdk4/6 and inactivates it; (11) RB1 binds to E2F inactivating it; (12) cyclinD1:Cdk4/6 phosphorylates RB1; (13) RB1-p binds to E2F; (14) cyclinE:Cdk2 phosphorylates RB1-p, preventing it from binding E2F; (15) free E2F enhances production of c-Myc; (16) free E2F enhances production of RB1; (17) free E2F drives the G1–S transition and proliferation. (b) Mathematial model wiring diagram. The species in red were measured by experiment and those in black were not. For p21, cyclinD1 and RB1, the total protein was measured. In addition, RB1-pp was measured.

In the mathematical model (figure 1b), we make some simplifying assumptions to reduce the number of species modelled so as to be more in line with the number of species measured. In particular, we do not model Cdk4/6 explicitly, but assume that all cyclinD1 not bound to p21 is bound to Cdk4/6 and active. Thus, to model the effect of palbociclib, which inactivates cyclinD1:Cdk4/6, we allow palbociclib to bind to cyclinD1 in the model and hold it inactive. Similarly, we do not explicitly model Cdk2. In addition, we do not model E2F, but assume that the level of RB1-pp reflects the transcriptional activity of E2F. While not modelling E2F may seem a step too far, the actual biological complexity is more significant than that shown in figure 1a. There are three pocket proteins (RB1, p107 and p130) that bind to the nine different E2F species (E2F1–9), and the binding efficiencies are governed by many phosphorylation sites on the pocket proteins [24]. Since our ultimate objective is predicting proliferation, we settled on associating the RB1-pp level (associated with phosphorylation on S612, which is associated with proliferation [25]) with proliferation, as this matches our experimental data.

Modelling oestrogen deprivation

Removing E2 from cultured cells that have been growing in medium containing E2 cannot be accomplished by simply changing to a medium containing no E2 [26]. MCF-7 cells growing in +E2 conditions have a much higher internal concentration of E2 than that of the medium due to non-specific binding of E2 in the cytoplasm as well as specific binding of E2 to various oestrogen receptors in the cell. When the medium is changed to a –E2 medium, E2 from the cells leaches into the new medium, establishing a new E2 level that can be significant for maintaining proliferation as MCF-7 cells can proliferate significantly in 10 pM of E2 [27]. In our 7-day experiments, the medium was changed to treated medium on day 0 and on day 3. For treatments involving E2 deprivation, we incorporated two additional model parameters, one for the E2 concentration in the medium after day 0 and one for the concentration after day 3, to account for the fact that the E2 concentration does not go to zero and that it decreases with additional changes in the medium.

Parameter estimation

The degradation rates for different proteins were assigned according to half-lives found in the literature, where k = −log(1/2)/Thalf-life. The other parameters were optimized to minimize the discrepancy between the model simulation and the experimental result. The cost function waswhere i indexes the state variables (proteins or cell number), j indexes the experimental time points, is the experimental measurement of the ith variable at time j, y(t, p) is the simulation result of the ith measurement at time j using parameter vector p, and σ is the standard deviation of the measured value based on three replications. The data and time points used for parameter estimation are listed in table 1. Parameter estimation was performed using MATLAB [28]. Initially, we tuned the model manually to get the gross responses correct. Then we relied on the default genetic algorithm function, ga, in MATLAB's global optimization toolbox and used the fminsearch function to refine the results of ga. We broke the parameters into groups, initially tuning the parameters related to ER, then tuning those related to all other proteins, and finally tuning all these parameters together. The three parameters associated with palbociclib were tuned separately after the initial model was tuned. The parameters associated with proliferation were tuned last, and finally all the parameters were tuned together to produce the final result.

Statistical analysis

The Mann–Whitney U-test was used for statistical comparisons [29]. The western blot protein data from the treatment conditions were first normalized to the control condition and then compared with 1 to test whether the treatments significantly decrease or increase the protein levels after time 0. Cell numbers after treatment were directly compared with the control numbers to test whether the treatments decreased cell proliferation.

Results

The proposed model structure can largely explain the experimental data

The training data for estimating the model parameters consisted of time-course measurements of the proteins in red in figure 1b for the –E2 and +E2 + ICI treatment conditions. These data are shown in figure 2, and it can be seen that the majority of measurements are statistically significant (asterisks), although the p21 measurements are quite noisy, and there is no real trend in the cyclinD1 data for the –E2 case. However, there are clear trends in most of the data, and it is critical that the model structure is capable of reproducing these trends. Figure 2a shows the experimental observations over 7 days (red) in response to –E2. The significant ER increase after treatment can be captured by the model because the half-life of E2:ER (oestrogen-bound ER) is about 3–4 h [30] and the half-life of unbound ER is about 4–5 h [30]. Therefore, the depletion of E2 stabilizes ER and causes the total ER level to increase. The jump in ER after 3 days is due to the medium exchange that decreases the concentration of E2 in the medium that is leached from the cells (see Methods). The small decrease in the cyclinD1 level can be captured by the model since E2:ER is a transcription factor for cyclinD1. By contrast, c-Myc not only decreases at the beginning, but decreases abruptly at 3d and continues to decrease over the remainder of the time course. The model captures the sharp decreases due to changes in the E2 concentration levels at 0d and 3d and the gradual decrease via RB1-pp, which drives the production of c-Myc. As the cells are gradually blocked in G1, RB1-pp decreases steadily, causing a decrease in c-Myc. The noisy p21 data do appear to increase; this is captured in the model as decreasing c-Myc releases its transcriptional inhibition of p21. We should note that there are likely to be many ways in which c-Myc influences RB1-pp, such as activation of cyclinE:Cdk2 kinase activity via Cdc25A [31], but, for simplicity, we have chosen to model c-Myc's effect only through the p21 pathway. The model can capture the reduction in RB1-pp, since this follows from the reductions in c-Myc and cyclinD1.

Figure 2.

Simulation and experimental results for protein-level changes under the different treatment conditions used for model calibration. The experimental measurements are shown in red (mean value ± s.e., n = 3) and the simulation results are shown in cyan (lowest cost simulation as a solid line and the interval containing the central 98% of the cohort simulations as a shaded region). Statistically significant changes from T = 0 are denoted by an asterisk. (a) Protein level changes after E2 deprivation over 7 days. (b) Protein-level changes after +E2 + ICI treatment over 7 days. ICI concentration is 500 nM and +E2 concentration is 10 nM. The simulation values for total cyclinD1 and total p21 are plotted to enable comparison with experiment.

At first glance, it seems odd that total RB1 should decrease much more than RB1-pp and that this decrease should be associated with the observed decrease in proliferation, as less total RB1 might be expected to create free E2F and promote the G1–S transition. This decrease in total RB1 has been observed by others in response to E2 deprivation [32]. In the model, free RB1 is always much more prevalent than either of the phosphorylated forms. The level of RB1-p is controlled by the cyclinD1 level, since the influence of free RB1 is saturated in the Hill function kinetics and has little effect on the RB1-p level. In turn, the level of RB1-pp is controlled by cyclinE and not influenced significantly by the RB1 level. The model captures the changes in total RB1 through the presumed transcriptional effect of RB1-pp on RB1.

Figure 2b shows the experimental observations in response to treatment with ICI. The model can capture the decrease in ER since ICI binds to ER, displacing E2 and enhancing ER degradation via the proteasome [33]. The reduced availability of the transcription factor ER:E2 due to ICI treatment causes changes in c-Myc, cyclinD1, p21, RB1-pp and total RB1 similar to those caused by reduced availability of ER:E2 due to E2 deprivation, and for the same reasons.

Model calibration and the simulation cohort

To calibrate the model, we used a genetic algorithm to search for parameter values minimizing the least-squares cost function of §2.7. To address the fact that the model may not be practically identifiable from the data, we identified 399 other parameter sets that provided a reasonable fit to the data to use as a simulation cohort. Any member of the cohort is a reasonable parameterization of the model and when we use the model for prediction we simulate all 400 parameter sets. The resulting spread in the predictions enables us to see how well the data used to calibrate the model constrain the prediction. The coefficients of variation of the parameters and the variations in the cost function for the simulation cohort are plotted in electronic supplementary material, figure S3.

The final results of calibrating the model to the training data are shown in figure 2. The results are quite reasonable, given the simplicity of our model. The spread of results from the simulation cohort reflects both the variation in the parameters for the chosen cohort as well as the sensitivity of the model to these parameters. The least spread is exhibited by cyclinD1 and RB1-pp. Although RB1-pp is not extremely sensitive to the parameters driving it (see sensitivity analysis in the electronic supplementary material), proliferation is very sensitive to RB1-pp level and forces its level to be heavily constrained. CyclinD1, on the other hand, is simply not very sensitive to the parameters driving it and the limited coefficients of variation of these parameters in the cohort result in a small spread.

Of course, there is always the question of whether the fit is reasonable because the model captures the essentials of the mechanism or simply because enough parameters were added to the model. In the following, we argue that our ability to add a new drug to the model without recalibration, and the ability of the model to predict drug combinations, shows that the model is capturing some of the essential mechanism.

Adding a new drug to the model requires limited new data

If the model is capturing the basic mechanism driving protein changes, adding a new drug to the model should require only calibrating a few new parameters associated with the new drug. Such calibration should not require measuring the response of all proteins in the model to the new drug, but only a few key proteins. We illustrate this by incorporating palbociclib, which has been used clinically in combination with endocrine therapy [5-7,10,34], into our model. Palbociclib inhibits Cdk4/6 kinase activity, reducing RB1 phosphorylation, primarily at S780/S795 [7,35], and ultimately leading to reduced hyperphosphorylation of RB1 and cell cycle arrest. To calibrate the binding and unbinding parameters associated with palbociclib in our model, we measured only two proteins, c-Myc and RB1-pp (S612, associated with proliferation but not directly affected by palbociclib [25]), which should experience strong downstream reactions to the drug. Figure 3 shows the results of calibrating the model to 1 μM of palbociclib added to the +E2 medium. As expected, the inhibition of cyclinD1:Cdk4/6 kinase activity by palbociclib decreases the RB1-pp level. This, in turn, leads to decreased transcription of c-Myc, causing the c-Myc level to decrease.

Figure 3.

Simulation and experimental results for protein-level changes in response to palbociclib. Simulations are shown in cyan and experimental measurements in red (mean value ± s.e., n = 3). Statistically significant changes from T = 0 are denoted by asterisks.

The proliferation results can be explained by the RB1-pp level

Proliferation is a primary endpoint for treatment decisions, so it is important for the model to be able to capture how proliferation changes with treatment. When considering treatments that produce G1 arrest, the key transition governing the execution of the cell cycle is the transition from the G1 phase to the S phase. A major determinant of the transition is the phosphorylation status of RB1 [36], but there are many other factors affecting the transition (see the discussion in the electronic supplementary material under Biological Justification of the Model).

Based on our experimental data, we assume that the rate of transition, hence of proliferation, in our model of an asynchronous population of cells is governed by the level of RB1-pp according toThe rate of cellular proliferation is modelled as proportional to the current number of cells subject to a carrying capacity constraint (logistic growth). The proportionality constant depends on the RB1-pp level via a Hill function. We expect proliferation to be near maximal at the pre-treatment RB1-pp level and, in our experimental data, proliferation essentially stops when the RB1-pp level reaches half of its normalized value, necessitating a relatively large Hill exponent of 6. The carrying capacity, assumed constant across all experiments, is used to account for the fact that in the control condition the cells approach confluence in 7 days; it has little effect on the other, slower growing, conditions.

Figure 4 shows that the model does a good job of matching the experimental proliferation results for the untreated (+E2), deprived (–E2), ICI-treated and palbociclib-treated cases on which the model was calibrated. This helps validate our claim that proliferation in these cases can be modelled using the RB1-pp (S612) level.

Figure 4.

Model simulations and experimental measurements of the normalized cell numbers under different treatment conditions. Experimental counts are in red (mean value ± s.e., n = 3) and simulations are in cyan. ICI concentration is 500 nM, palbociclib concentration is 1 µM and E2 concentration is 10 nM. Cell numbers are normalized to T = 0 value. Statistically significant changes from the +E2 control are denoted by asterisks.

The model can predict the effect of combination treatments

To partially validate the model, we test its ability to predict the effect of combination therapies on protein levels and cell proliferation. Figure 5 shows how the model simulations of combining E2 deprivation and ICI compare to the experimental results. With the exception of cyclinD1, the predictions match the experimental results for the six measured protein species well. By reducing the supply of E2 and the ER level simultaneously, larger changes in the protein levels are observed than with either monotherapy, as would be expected.

Figure 5.

Prediction simulations for protein-level changes in response to combination E2 deprivation +ICI therapy. ICI concentration is 500 nM. Experimental results are in red (mean value ± s.e., n = 3), simulations are in blue and monotherapy experimental results are in cyan. Statistically significant changes from T = 0 are denoted by asterisks.

The combination of E2 deprivation and palbociclib attacks a key G1 kinase, cyclinD1:Cdk4/6, by both reducing the level of cyclinD1 and inactivating Cdk4/6. As expected, the combination reduces the level of c-Myc and RB1-pp to a greater extent than either monotherapy. The predictions of the model are compared with experimental protein measurements in figure 6a. The agreement between measured and predicted c-Myc is good, but the simulation misses the low level of RB1-pp at the final two time points. Since proliferation in the model essentially stops when the normalized RB1-pp level falls below 0.5, it is still possible to capture the proliferation well while somewhat overestimating the RB1-pp level.

Figure 6.

Comparison of model predictions with experiment for combination therapies. Experimental data (red), simulations (blue) and monotherapy experimental results (cyan). (a) Prediction of c-Myc and RB1-pp levels in response to –E2 + palbociclib (1 μM) combination treatment. (b) Prediction of proliferation in response to –E2 + ICI and –E2 + palbociclib (1 μM) combination treatments. Asterisks denote statistically significant changes from T = 0 for proteins, and from the +E2 control for proliferation.

Figure 6b shows that the model also predicts the decreased proliferation in response to the combination therapies quite well. This follows from the fact that RB1-pp decreases more in response to the combinations than to either monotherapy.

Local sensitivity analysis of protein levels and proliferation

A local sensitivity analysis is used to check if the output of our model is very sensitive to certain parameters. The sensitivity coefficient for a given model output and a given parameter is the per cent change in the model output divided by the per cent change in the parameter value. To calculate this coefficient, we changed the parameter by ±5% from its best-fit value, so the sensitivity coefficient is not unduly local. In figure 7, we plot the sensitivity coefficients for the most sensitive outputs in our model, c-Myc, total RB1 and cell number, to each model parameter at the day 7 time point. The sensitivity coefficients of the other proteins are shown in the electronic supplementary material and are of the order of 4 or less, indicating relatively low sensitivity.

Figure 7.

Local sensitivity analysis of the effect on c-Myc level, total RB1 level and cell number of each model parameter. Each of the parameters is changed by ±5% and the sensitivity at day 7 for each treatment is plotted.

The most significant sensitivity for c-Myc occurs for the ICI treatment case and involves parameters related to the basal translation of cyclinE and p21 and the phosphorylation of RB1-p. This is not surprising as these parameters regulate the level of RB1-pp, which helps drive the c-Myc level in our model. Total RB1 is most sensitive for the –E2 + ICI case, and the parameters producing the greatest sensitivity involve basal transcription of cyclinE and c-Myc, the downregulation of p21 by c-Myc, phosphorylation of RB1-p by cyclinE and the dephosphorylation of RB1-pp. These parameters all converge on the level of RB1-pp.

Cell number is most sensitive in the ICI and palbociclib treatment cases and its greatest sensitivity is to parameters involving basal transcription of cyclinE and its phosphorylation of RB1-p. Again, these parameters drive the level of RB1-pp, which in turn drives proliferation. The fact that cell proliferation is essentially modelled as exponential accounts for the large sensitivity.

The model can be used to explore the effect of sequential therapies

While the initial response of patients to a mono- or combination therapy is often promising, in many cases, resistance to continuous therapy eventually arises (e.g. palbociclib [10,11], letrozole [8,12] and ICI [37]). The resistance mechanisms are varied, and in many cases not well understood, but can include protein overexpression, epigenetic changes, gene mutation and gene amplification or deletion [8-12,38]. Our idea to address some of these mechanisms is straightforward: alternate various therapies, perhaps including drug holidays, in a fixed cycle to significantly suppress proliferation of tumour cells, but not target one particular mechanism non-stop to avoid selecting for resistance mechanisms. The hope is that such an approach can hold off, or at least delay, the onset of resistance.

Optimizing sequential therapies under various restrictions is likely to be impractical experimentally because of the large number of variations in dosing and timing of multiple therapies that must be considered. On the other hand, a reliable mathematical model may be able to sort through the huge number of possibilities to find particularly promising therapy protocols that can be tested experimentally.

To illustrate the basic idea with our current model, we consider resistance to palbociclib, which is frequently associated with the overexpression of Cdk6, which increases over time with continuous palbociclib therapy [39,40]. Presumably, the increasing Cdk6 eventually binds enough palbociclib to free up cyclinD1:Cdk4/6 to phosphorylate RB1 and move the cell toward proliferation. In our model, we did not include Cdk4/6 explicitly but rather assumed that all free cyclinD1 was complexed with Cdk4/6 and, therefore, active. Palbociclib inhibits this activity in the model by binding to cyclinD1 and holding it inactive. We mimic the resistance effect of increased Cdk6 sequestering palbociclib by slowly increasing the level of cyclinD1 transcription under palbociclib therapy and allowing it to slowly decrease when palbociclib is withdrawn. So palbociclib is sequestered by cyclinD1, rather than Cdk6, in the model. Modelling details are in the electronic supplementary material, equations 40 and 41. Figure 8a shows how adding this resistance mechanism changes the predicted results of monotherapy, with resistance emerging after 30 days.

Figure 8.

Model simulations of possible alternating therapies compared with a continuous monotherapy with an added resistance mechanism. (a) Continuous palbociclib therapy with and without the resistance mechanism are compared. (b) Continuous palbociclib therapy is compared with alternation of E2 deprivation and palbociclib therapies. (c) Continuous palbociclib therapy is compared with alternation of ICI and palbociclib therapies. Doses: E2 (10 nM), ICI (500 nM), palbociclib (1 μM). The simulations use the lowest cost parameter set, and all simulations in (b) and (c) include the resistance mechanism.

Figure 8b compares a continuous therapy of palbociclib with a therapy that alternates E2 deprivation with palbociclib for one-week periods in a repeating cycle. While the cyclinD1 and RB1-pp levels slowly increase under the alternating therapy, the onset of resistance is effectively delayed. Figure 8c shows a similar effect when treatment with ICI and palbociclib are alternated with one-week periods.

Discussion

In this study, we created a relatively simple mathematical model to capture the effects of ER signalling on key cell cycle proteins governing progression through the G1–S checkpoint in breast cancer cells. The model was calibrated using experimental monotherapy data and was capable of predicting the experimental data from combination therapies. We further showed that it was relatively easy to add in the effect of a new drug that impinged on this checkpoint without having to measure the response of all the proteins in the model to the new drug or to recalibrate the previously determined parameters of the model. Finally, we linked the level of hyperphosphorylated RB1 (S612) to the proliferation rate and showed that this linkage could capture the proliferation observed in our experimental data.

Clearly, our model is a highly simplified version of the actual system. We assume that the behaviour of a population of cells can be adequately modelled by allowing the average levels of species in the population to interact according to a simplified version of the mechanism existing in individual cells. While there are theoretical objections to such a model (the interaction rate between A and B in a population is not necessarily proportional to the average concentration of A multiplied by the average concentration of B), this approach has produced useful models in other contexts [13]. One consequence of this approach is that our model, unlike the model of a single cell, will not exhibit a bistable switch governing the G1–S transition [17], as a population will not completely stop cycling even though individual cells will. Probably, the fact we are trying to model the behaviour of a large population of asynchronous cells allows the model to be simpler than if we were trying to model the intricacies of an individual cell.

Our ultimate goal is to use the model to derive alternating therapies that may stave off the development of resistance often associated with continuously applied therapies [41,42]. We illustrated some possibilities along this line, but clearly additional experimental confirmation is required. Confirmation that the model can accurately predict proliferation over longer time scales and in response to switching therapies will be the subject of future work, as will confirmation that cells subjected to alternating therapies are truly less resistant than cells subjected to continuous therapy.

We should also note that the treatments we are currently considering all impinge on the G1–S transition. Thus, a resistance mechanism that completely eradicates the G1–S checkpoint, such as loss of RB1, is likely to be unaffected by our approach. Drugs that arrest the cell cycle at a different point, for example G2, or that target pathways resulting in cell death will be necessary to address such cases. Drugs that induce significant cell death will also be necessary if we want to provide a drug holiday. While the holiday may help prevent the development of resistance, it will also allow significant proliferation and holding down the net proliferation will require a drug that can kill cells [43]. Finally, we realize that successful results in vitro are potentially a long way off from success in animal models or human tumours, but it is the place to start for mechanism-based modelling.