## Abstract

Endocrine therapy, targeting the oestrogen receptor pathway, is the most common treatment for oestrogen receptor-positive breast cancers. Unfortunately, these tumours frequently develop resistance to endocrine therapies. Among the strategies to treat resistant tumours are sequential treatment (in which second-line drugs are used to gain additional responses) and intermittent treatment (in which a ‘drug holiday’ is imposed between treatments). To gain a more rigorous understanding of the mechanisms underlying these strategies, we present a mathematical model that captures the transitions among three different, experimentally observed, oestrogen-sensitivity phenotypes in breast cancer (sensitive, hypersensitive and independent). To provide a global view of the transitions between these phenotypes, we compute the potential landscape associated with the model. We show how this oestrogen response landscape can be reshaped by population selection, which is a crucial force in promoting acquired resistance. Techniques from statistical physics are used to create a population-level state-transition model from the cellular-level model. We then illustrate how this population-level model can be used to analyse and optimize sequential and intermittent oestrogen-deprivation protocols for breast cancer. The approach used in this study is general and can also be applied to investigate treatment strategies for other types of cancer.

## 1. Introduction

Breast cancer remains one of the most common cancers in women. Approximately 70% of breast cancers express oestrogen receptor-α (ERα), which is a major driver of survival and proliferation in breast cancer cells [1]. Endocrine therapies, targeting the ERα pathway, have become routine and effective treatments for breast cancer [2]. Three types of endocrine therapies are commonly used in breast cancer: (i) aromatase inhibitors (AIs) to deprive cells of 17β-oestradiol (E2, the primary oestrogen in breast tumours), by reducing E2 biosynthesis, so as to remove the ligand of ERα; (ii) selective oestrogen receptor modulators (SERMs, e.g. tamoxifen), to inhibit ERα binding with E2; and (iii) selective oestrogen receptor downregulators (SERDs, e.g. Faslodex), to inhibit E2 binding and target the ERα protein for degradation. Despite their success in reducing the annual death rate from breast cancer by one-third over the past decades, resistance to endocrine therapies often develops, creating a fundamental problem in breast cancer treatment [2].

The development of resistance to endocrine therapy has often been thought of as a progressive, step-wise phenomenon whereby cancer cells gradually shift from an E2-dependent, responsive phenotype to a less-responsive phenotype, and ultimately to an E2-independent phenotype [3,4]. This progression implies that endocrine treatments should be sequenced to take advantage of the cancer cell's oestrogen-responsiveness state. Pre-clinical and clinical observations suggest that breast cancer with acquired resistance may respond to the second- and third-line endocrine treatments (figure 1*a*). For example, a proportion of breast cancers that have developed resistance to tamoxifen can still be inhibited by AIs, and some AI-resistant breast cancers can still benefit from Faslodex [3]. Nonetheless, cross-resistance among these therapies occurs, and both the frequency and duration of responses to the second- and third-line therapies are often less than those achieved with a first line intervention. Various combinations of consecutive anti-oestrogen treatments are still under clinical investigation [6]. In addition, giving breast cancer patients a ‘drug holiday’ has been suggested to deal with acquired endocrine resistance [7–9] (figure 1*a*). The duration of treatment and the interval between treatments are often critical parameters for the effectiveness of any intermittent treatment regimen. However, a rigorous understanding of how breast cancer patients can benefit from intermittent endocrine treatment is still lacking, preventing us from optimizing treatment strategies. Optimization of sequential treatment faces a similar problem [10].

In this initial study of devising a time-varying therapy to prevent the onset of resistance, we consider only the effect of E2-deprivation therapy for simplicity. A mathematical model is built to account for the different E2-responsiveness phenotypes observed in breast cancer cells and used to explore the transitions from endocrine responsiveness to endocrine resistance under sequential and intermittent treatment.

It is widely reported that there are three phenotypes (states) of breast cancer cells related to E2 sensitivity [3,5] (figure 1*b*). Both clinical observations [11] and cell line studies [12] suggest that breast cancer cells, which are E2 sensitive initially (figure 1*b*, left), can adapt to E2 deprivation by entering an E2-hypersensitive state, for which a non-physiologically low concentration of E2 is sufficient to support cell proliferation (figure 1*b*, middle). Prolonged and profound E2 deprivation can lead to an E2-independent state [3,5] (figure 1*b*, right).

Notably, not all cells make the state transitions necessary for proliferation. Rather, the E2 sensitivity progression described earlier most probably represents somatic evolution, whereby cells that switch to the hypersensitive or independent states outgrow their disadvantaged peers. The diversity of cell phenotypes that is necessary for somatic evolution can be caused by genomic or non-genomic heterogeneity [13]. In this paper, we consider only non-genomic heterogeneity, whose occurrence is supported by experiments that demonstrate the reversibility of E2-hypersensitive and -independent phenotypes in breast cancer cells during two years of E2 manipulation *in vitro* [14]. We do not address the issue of genomic heterogeneity due to mutation, which probably occurs over a much longer timescale.

In physics, an intuitively appealing way to visualize the dynamics of a system is to view the system as a particle moving on a surface defined by a ‘potential energy’. In two dimensions, this surface can be plotted as a landscape where the particle moves downhill and falls into basins. When the landscape contains multiple basins of attraction, there are multiple stable states where the particle can remain indefinitely. If the system is subjected to random external forces, it is possible for the system to transition from one basin of attraction to another, with important implications for phenotypic plasticity. Inspired by this idea C. Waddington proposed describing cell fate decisions by the dynamics of a control system on an ‘epigenetic landscape’ [15]. The basins of attraction were (in Waddington's view) stable cell fates, and the paths that cells take across the landscape corresponded to the differentiation events of cells in a developing organism. In recent years, the ability to compute a potential landscape from an underlying mathematical model has been used to visualize a variety of cell fate decisions [16–22]. While the imagery of Waddington's landscape has been used to discuss the fate of cancer cells during the development of resistance to therapy [18,21], little has been done to quantitatively map this landscape from a mathematical model of the underlying signalling network of cancer cells.

Based on the mathematical, cell-level model that accounts for the three E2 sensitivity states observed in breast cancer cells, we compute the potential landscape, which we refer to as the oestrogen response landscape. This landscape provides a visual picture of how breast cancer cells change phenotypes in response to variations in protein expression (noise) and varying environmental conditions (E2 withdrawal therapy). We also evaluate the effects of population selection in reshaping the oestrogen response landscape. In addition, we apply a multiscale modelling strategy to estimate the transition probabilities in a simple population-level, state-transition model, based on the minimum action path (MAP) between basins of attraction on the oestrogen response landscape [19,23]. By using this state-transition model, we investigate the population dynamics of breast cancer cells under both sequential and intermittent treatment. The model shows that a proper treatment sequence is important to achieve maximum response and to avoid cross-resistance. The quantitative model also allows us to search for an optimal protocol for intermittent treatment to overcome acquired endocrine resistance. Thus, our study provides a methodology for exploring breast cancer treatment strategies grounded in cellular-level, biologically based, mathematical models.

## 2. Material and methods

### 2.1. Three oestrogen receptor-α activation modes and crosstalk with growth factor receptors

The ERα pathway in breast cancer cells has at least three distinct branches [2,24–26]: (i) the E2 branch, where ERα bound to E2 becomes an active transcription factor of genes supporting cell survival and proliferation; (ii) the phosphorylation branch, where ERα becomes an active transcription factor subsequent to phosphorylation at multiple sites (including Ser118, 104, 106 and 167) in its AF1 domain and (iii) the membrane branch, where membrane-bound ERα mediates the activation of the MAPK and PI3K/AKT pathways (the ‘non-genomic’ role of ERα, since it is not acting as a transcription factor in this branch) [2]. These three branches of the ERα pathway, encompassing both the genomic and non-genomic roles of ERα, crosstalk with each other and with the growth factor receptor (GFR) signalling network (figure 1*c*) [27,28]. We consider the following four crosstalk mechanisms: (i) E2-liganded ERα (E2ER in figure 1*c*) inhibits expression of both GFR [29,30] and membrane ERα [31] (ERM in figure 1*c*); (ii) both GFR and ERM activate downstream kinases that phosphorylate ERα [2,24] (ERP in figure 1*c*); (iii) ERM has positive interactions with GFR signalling at different levels [2,24]; (iv) both ERM and GFR are involved in auto-activation loops at the signalling and/or epigenetic levels [5,31–33]. Experimental data indicate that the different E2 sensitivities in breast cancer cells are determined by the balance of the three ERα signalling branches (E2ER, ERP and ERM) and GFR signalling [3,5]. In particular, E2-sensitive cells depend on E2ER for proliferation, whereas E2-hypersensitive cells depend upon active ERM, which contributes to a high level of ERP. In E2-independent cells, proliferation depends on high GFR activity, which also produces a high level of ERP. These experimental observations are summarized in figure 1*d*.

We create a phenomenological model that of necessity does not attempt to incorporate complex mechanistic details, many of which are inadequately known, but rather to reproduce the experimental observations summarized in figure 1*d* using the simplified signalling network of figure 1*c*. Since ERP activity does not feedback on production or destruction of other species in the model, and its state is determined by the state of ERM and GFR, we do not include ERP in our dynamical model. In addition, we assume the binding of E2 to ERα is sufficiently fast that the amount of E2ER can be modelled as an algebraic function of the E2 level. Hence, we model only the dynamics of ERM and GFR, incorporating the influences shown in figure 1*c*. Proliferation is modelled by using figure 1*d* and the activity levels of ERM and GFR to determine a cell's E2-sensitivity state. The cell is then assigned a proliferation index (PI), as a function of E2, based on this sensitivity as in figure 1*b* (see also figure 2*a,b*). This procedure reproduces the experimental observations that E2-sensitive cells grow in high levels of E2, but die in low or trace levels, whereas E2-hypersensitive cells grow in high or low E2 concentrations, but die in trace levels, and E2-independent cells grow under all E2 conditions. The numerical rates of growth and death, however, are not based on experiment.

### 2.2. Model implementation

Our mathematical model of the influence diagram in figure 1*c* uses stochastic differential equations (SDEs) to capture both the deterministic effects of species' interactions in a signalling network and the random effects of intrinsic and extrinsic sources of noise [34–37]. Let X* _{i}* be the name of the

*i*th species in our interaction diagram (E2ER, ERM and GFR). Letting [X

*] = activity of X*

_{i}*and*

_{i}*x*= log

_{i}_{10}[X

*], we write the SDEs for our model as 2.1where 2.2We assume that all concentrations in our model have a dynamical range of about 10-fold, so, on a log*

_{i}_{10}scale, the

*x*'s vary between 0 and 1. The parameters

_{i}*γ*determine the rates at which the variables

_{i}*x*approach their time-varying ‘target’ values,

_{i}*x*=

_{i}*H*(

*W*). Note that

_{i}*H*(

*W*) is a sigmoidal function that varies between 0 and 1, with a value of 1/2 when

_{i}*W*= 0.

_{i}*W*is the net activation or inhibition on species X

_{i}*, and the leading term,*

_{i}*ω*, determines whether X

_{i}*is activated or inhibited when there are no regulatory signals impinging on X*

_{i}*from any species in the network. In the summation,*

_{i}*ω*is the strength of the interaction from X

_{ij}*to X*

_{j}*(*

_{i}*ω*> 0 for activation and <0 for inhibition).

_{ij}*ζ*(

_{i}*t*) is a temporally uncorrelated, statistically independent Gaussian white noise process with 2.3where and

*D*characterizes the magnitude of the random fluctuations on each species.

_{i}The details of our model are provided in the electronic supplementary tables. Electronic supplementary material, table S1 defines the variables of our model, table S2 lists the precise equations of the model and table S3 records the parameter values that we used for our simulations. In deriving the model equations, we applied a quasi-equilibrium approximation to ERE2*,* because the binding of E2–ER occurs on a much faster timescale than variations in ERM and GFR, which also may have transcriptional, possibly epigenetic, regulation. Hence, only ERM and GFR are described by SDEs. We also assumed that ERM and GFR are subjected to similar random fluctuations, so *D*_{1} = *D*_{2} = *D* in equation (2.3). We used Matlab (v. 7.9.0) to perform simulations.

Our modelling is done at a phenomenological level; consequently, none of our parameters are directly related to measurable physical rate constants of the system. Phenomenological modelling is both reasonable and necessary, because the rate constants for most of the interactions in the model are unknown. The parameters are assigned to match the phenotypic properties of the ERα–GFR control system. Initially, the parameters for the deterministic part of the model were chosen to ensure that the system has attractors corresponding to the E2-sensitive, -hypersensitive and -independent states (figure 2*a*). For the stochastic simulations, we chose the noise parameter, *D* = 0.005, to ensure that the system's transitions occurred in a manner similar to experimental observations. To model selection pressures, we chose the shape of the PI curves for cells in each state (figure 2*b*) to capture qualitatively the experimentally observed growth properties of E2-sensitive, -hypersensitive and -independent cells.

### 2.3. Mapping the oestrogen response landscape

A dynamical system whose state moves in the direction of the negative gradient of a potential function can be conveniently visualized by plotting the potential energy landscape and viewing the system state as a particle that moves downhill on this landscape, coming to rest at a local minimum. While open biological systems, such as the one we consider here, are not governed by the gradient of a potential function, the ideas in [17–20,38] can be used to compute a potential landscape useful for visualization. The basins of this landscape are regions of high probability for finding the system, and the hills are regions that the system rarely visits.

We write the governing equation (2.1) of our breast cancer model in vector form as
2.4where , with similar definitions of **F** and **ζ**. Let *P*(**x**,*t*) be the probability distribution function for state **x** at time *t*, i.e. *P*(**x**,*t*)d*x*_{1}d*x*_{2} = probability of finding the system in the small area (*x*_{1}, *x*_{1} + d*x*_{1}) × (*x*_{2}, *x*_{2} + d*x*_{2}) at time *t*. The equation governing *P*(**x**,*t*) is
2.5where and
2.6is the probability flux for our system.

We wait long enough for the probability density function, and hence the flux, to reach a steady state, and then solve equation (2.6) for **F**(**x**) in terms of *P*_{ss}(**x**) and **J**_{ss}(**x**)
2.7where *U*(**x**) = −ln(*P*_{ss}(**x**)). Thus, **F** can be expressed as the gradient of a potential function *U*(**x**) plus a flux field. Although the system's dynamics is not governed solely by the gradient of *U*(**x**) but includes as well a ‘spiralling’ motion due to the flux field [15], the generalized potential *U* = −ln(*P*_{ss}) still provides a useful visualization of the dynamics of the system. We refer to *U*(**x**) as the oestrogen response landscape in this study. Because of how *U*(**x**) is related to *P*_{ss}(**x**), basins of *U*(**x**) correspond to high probability states and hills to low probability states. Details for solving the diffusion equations are provided in the electronic supplementary material, Text S1.

### 2.4. Cell state transition model

In the small noise case, the system spends most of its time fluctuating around one of its *m* steady states (*m* = 4 in this study) with occasional transitions between states. Hence, it is inefficient to track the states of cells in a large population by simulating each cell individually. To study state switching within a large population of cells, we transform the cellular-level, multiple steady-state model into a Markov state-transition model by using the cellular model to compute the Markov state transition rates. Hence, we divided the two-dimensional state space of the cellular-level model into four blocks, representing the four discrete phenotypes, which we refer to as ‘regions’ R_{1}, R_{2}, R_{3} and R_{4} in figure 2*a*. Similar to [39,40], we describe a population of cells by a vector, **v**, of length *m* = 4. The *i*th entry of the vector is the number of cells in region R* _{i}*, and the total number of cells in the population at time

*t*is . The population vector

**v**evolves according to the equation where

**A**is an

*m*×

*m*transition rate matrix. The off-diagonal elements of

**A**,

*A*=

_{ij}*k*, are the transition rates from region

_{ij}*j*to region

*i*, and the diagonal elements are , where

*u*is the proliferation rate of cells in region

_{i}*i*. Thus, starting from a vector reflecting the initial populations in R

_{1}, … , R

_{4}at time 0, we can easily compute the temporal evolution of the population (for details, see the electronic supplementary material, Text S2). The average PI of the population at each time step is then calculated by . PI is the specific reproduction rate of the population: PI > 0 signifies an expanding population and PI < 0 signifies a diminishing population.

### 2.5. Minimum action path

To estimate the transition rates *k _{ij}* among the different phenotypic regions of the dynamical system, which we need for the population-level model, we followed the approach in [19,23,41] based on applying the Wentzell–Freidlin theory [42] to equation (2.4) for the ‘small noise’ case. The key idea of this theory is that the most probable transition path from R

*at time 0 to R*

_{j}*at time*

_{i}*T*, , minimizes the action functional (equation (2.8)) over all possible paths 2.8

This optimal path is referred to as the MAP. We are interested in MAPs between fixed-point attractors (stable steady states) in regions, R_{1}, … , R_{4}. The MAPs connecting each pair of attractors were computed numerically by applying the minimum action method used in [23,41]. For *T* large enough to allow a minimum action transition from the attractor in R* _{j}* to the attractor in R

_{i}, becomes independent of

*T*, which is the goal. In our case,

*T*= 40 is a reasonable choice: large enough but not too large (see the electronic supplementary material, figure S1).

The MAP is useful for our purposes because we find that *k _{ij}*, the transition rate from state R

*to state R*

_{j}*, is related to by the empirical equation 2.9where the parameters*

_{i}*α*,

*β*and

*n*are estimated by fitting equation (2.9) to the results of direct Monte Carlo simulations of the dynamic model, where the Langevin equation (2.4) is integrated by the explicit method 2.10Here the

*η*are independent, normal, random deviates with mean 0 and variance 1.

_{i}## 3. Results

### 3.1. The oestrogen response landscape

Both clinical observations and cell line studies suggest that breast cancer cells can shift between three different E2 sensitivity states [3,5] (figure 1*b*): E2-sensitive, -hypersensitive and -independent. In E2-sensitive cells, E2-liganded ERα (E2ER) is the major driver of pro-survival and pro-proliferation signalling, whereas ERα in the cell membrane (ERM), phosphorylated ERα (ERP) and GFRs are kept at low levels of activity. As summarized in [3,5], E2-hypersensitive cells are closely associated with elevated ERα in the cell membrane (ERM), where ERα plays a non-genomic role to mediate the activation of the MAPK and PI3K/AKT pathways. GFRs, such as IGF1R and HER2, are not necessarily overexpressed in E2-hypersensitive cells. A possible explanation for hypersensitivity is that ERα, phosphorylated by kinases in the upregulated pathway, may have increased affinity for E2. In E2-independent cells, however, increased GFR and possibly other signalling pathways play pivotal roles, probably by activating MAPK and PI3K/AKT that further phosphorylate ERα (ERP), driving E2-independent growth. These observations are summarized in figure 1*d*. The activity levels of ERM and GFR determine the E2-sensitivity state of our model. Thus, of the four regions in figure 2*a*, R_{1} corresponds to E2-senstivity (GFR low, ERM low), R_{2} corresponds to E2-hypersensitivity (GFR low, ERM high), and R_{3} and R_{4} correspond to E2-independency (GFR high, ERM high or low).

As discussed in Material and methods, the parameters of our model were chosen to reproduce the qualitative experimental data summarized in figure 1*d* [3,5]. The timescale of the model is not explicitly specified, but based on experimental observations we expect it to be on the order of weeks or months. The oestrogen response landscape computed from our model is plotted in figure 2*c* for three different E2 conditions ([E2] = 10^{2}, 10^{0} and 10^{−2}) corresponding to HIGH, LOW and TRACE amounts of E2. Varying the level of E2 reshapes the landscape and can induce cells to transition from one state to another. When E2 level is HIGH, the potential landscape has only one dominant basin in R_{1}, where both ERM and GFR are low (figure 2*c*, left), indicating an E2-sensitive state. Most cells will stay in this region. Reducing E2 to LOW level (figure 2*c*, middle) reshapes the landscape, deepening the basin in R_{2} and causing some cells to shift from the E2-sensitive state to the E2-hypersensitive state. This transition allows the oestrogen receptor to become activated at much lower E2 concentrations, essentially inducing a left-shift in the proliferation dose–response curve. Further decreasing E2 to TRACE amounts creates a dominant basin in R_{4}, the E2-independent state (figure 2*c*, right), and the E2-sensitive basin in R_{1} becomes very shallow. Both E2-hypersensitive and -independent cells will coexist in this condition, with a preference for the independent state.

In addition to cell state transitions (cellular adaptation) in response to a changing E2 environment, breast cancer cells are also under selection pressure due to the different proliferation rates that apply to each cell state in a given E2 environment. To model this selection pressure, the state space of the model was divided into four regions (R_{1}, R_{2}, R_{3} and R_{4} as before), with each region assigned a specific PI as a function of [E2] (figure 2*b*). Allowing cells in each region to proliferate at their assigned rates given the specified E2 concentration results in the landscapes in figure 2*c* (see the electronic supplementary material, Text S1 for calculation details). Selection pressure deepens some basins and destroys others (figure 2*d*). In particular, for LOW values of E2, selection pressure causes the dominant basin to shift from the E2-sensitive state (R_{1}) to the E2-hypersensitive state (R_{2}).

Thus, our model reproduces the different cell states observed experimentally in response to E2 manipulations. The oestrogen response landscape we have computed here provides a useful visualization of the global behaviours of the breast cancer system with three different E2 sensitivity states.

### 3.2. A state-transition model

To move from one basin to another, cell state transitions require the system to cross a barrier. The transitions are driven by intrinsic and extrinsic noise, with lower barriers corresponding to more probable transitions. The most probable routes for transitions among the four regions on the oestrogen response landscape, for a specified value of E2, were computed by finding the paths that minimize the action, a measure of the difficulty of moving from one region to another (see Material and methods), over all possible paths that start and end in the desired regions. An example of MAPs on the oestrogen response landscape is shown in figure 3*a*. Interestingly, the route taken from one basin to another is not the same as the route of the reverse transition, as seen by the non-overlapping yellow and red curves in figure 3*a*. This divergence of MAPs is a consequence of the non-gradient component of the flow field (equation (2.7)) [16,17,43].

The action along the MAP between two steady states in figure 3*a* can be used to estimate the transition probability between the two regions [19]. Thus, we again divide the state space into four regions (R_{1}–R_{4} in figure 3*b*) corresponding to the four cancer cell phenotypes. Each region has at most one attractor. For various E2 levels, the transition rates between pairs of regions were computed by Monte Carlo simulations of the stochastic model (electronic supplementary material, figure S2). In figure 3*c*, we plot these transition rates as functions of the action of the MAP between the corresponding fixed points. As expected, the points all lie close to a single curve described by equation (2.9), with *α* = 80.48, *β* = 0.07197 and *n* = 0.9805. Thus, for different values of [E2] and any pair of fixed points we can compute the action, and then use equation (2.9) to estimate the transition probability, without recourse to additional Monte Carlo simulations. The comparisons between transition rates estimated from actions and from direct Monte Carlo simulations at different E2 levels are shown in figure 3*d*,*e*. Because the transition rates out of R_{3} are very much greater than the transition rates into R_{3} (, ), the third block is rarely populated and can be safely neglected. The transition rates between R_{1} and R_{4} are much smaller than their transition rates with R_{2} (); hence, direct transitions between R_{1} and R_{4} can also be safely neglected.

Given the estimated transition rates, we built a deterministic model of phenotype transitions and simulated the growth of cancer cell populations under different E2 conditions. Figure 4*a* shows the steady-state proportion of each phenotype as a function of E2 level for the case of no selection pressure (PI = 0 for each phenotype). This distribution of phenotypes is consistent with the landscape plot in figure 2*c*. Note that when E2 is LOW ([E2] = 1), although there is a small population of E2-hypersensitive cells (green curve), the majority of cells are still in the E2-sensitive state (blue curve). However, adding the effect of population selection significantly changes the outcome (figure 4*b*). As in figure 2*d*, the selective advantage of E2-hypersensitive cells in the LOW E2 environment substantially increases the proportion of E2-hypersensitive cells. The proportion of E2-independent cells under TRACE E2 conditions ([E2] = 10^{−2}) is greatly enriched (figure 4*b*, right). Corresponding steady-state probability distributions, *P*_{ss}(*x*_{1},*x*_{2}), in the state space of the cellular-level model are shown in figure 4*c*,*d*. Simulating the phenotype transition model at the population level provides a substantial reduction in computation compared with Monte Carlo simulations of the cellular-level model. Consequently, we can easily study the consequences of both sequential and intermittent treatments of breast cancer under a variety of conditions.

### 3.3. Sequential treatment of breast cancer

The effects of sequential E2 deprivation on breast cancer cells have been reported in both clinical and pre-clinical settings [3,5,11]. For example, in pre-menopausal women, E2-dependent breast cancers often regress after surgical removal of the ovaries, which lowers circulating E2 from 200 to 15 pg ml^{−1}. After an initial stage of regression, breast cancer often regrows. AIs, now given as second-line therapies, can induce additional tumor regression by further decreasing E2 to less than 5 pg ml^{−1} [11].

Creation of an additional ‘response window’ using a second-line endocrine therapy can be captured by our state transition model (figure 5*a*). To illustrate this, we assume a primary therapy that can reduce the E2 level from HIGH to LOW, and a secondary therapy that can reduce the level from HIGH to TRACE. Decreasing E2 in a step-wise fashion (figure 5*a*, top), by first applying the primary therapy and then changing to the secondary therapy, results in sequential transitions of the cell population from the E2-sensitive to the E2-hypersensitive to the E2-independent state (figure 5*a*, middle). At each transition, cells in the previously favoured state begin to die off while some make the transition to the newly favoured state and begin to grow. Plotting the average population PI versus time clearly shows the two response windows of negative average growth (yellow phases, figure 5*a*, bottom).

Notably, if the E2 steps are applied in the reverse order (figure 5*b*), corresponding to applying the secondary therapy first, the simulation shows no second response window. This might be viewed as an example of the phenomena of cross-resistance to endocrine therapies in breast cancer. In figure 5*b*, the green phase represents the cross-resistance of cells to endocrine treatment associated with LOW E2 when the cells have been pretreated with (and gained resistance to) endocrine treatment associated with TRACE levels of E2.

### 3.4. Intermittent treatment of breast cancer

Intermittent E2 treatment to reverse acquired resistance and provide a second response window is simulated in figure 5*c*. Shifting between two different endocrine therapies can also produce a similar result (figure 5*d*). These simulation results provide support for both pre-clinical and clinical observations that intermittent treatment may be better than continuous treatment [7–9].

The state transition model can serve as a quantitative tool to find the optimum values for treatment duration, *T*_{treat}, and treatment breaks, *T*_{break}. Figure 6*a* shows the results for *T*_{treat} = 3000 and *T*_{break} = 2000, where cell population growth is reduced—but not stopped—by intermittent treatment. By contrast, intermittent treatment with *T*_{treat} = 1000 and *T*_{break} = 2000 (figure 6*b*) causes the total cell population to stop growing, thereby controlling—but not eliminating—the cancer. Comparing these two cases, we note that shorter treatment durations sometimes may be more beneficial than longer treatments.

Using the model to explore the effects of the two key treatment parameters on the expansion of the cancer cell population, we plot in figure 6*c* the total number of cells after treatment as a function of *T*_{treat} and *T*_{break}. Since the simulation starts with 10^{3} cells, any parameter values within the 10^{3}-contour (superimposed on figure 6*c*) will cause cancer progression to be arrested or reversed. These results point out the potential utility of a population-level, state-transition model grounded on a cellular-level, mechanistic model for optimizing cancer treatment protocols.

## 4. Discussion

As a heterogeneous population of cells, cancer can respond by adaptation to many changing stressors and stimuli in its environment. This plasticity can lead to acquired resistance and drug failure, a fundamental problem in cancer treatment. One widely accepted explanation for cancer plasticity is that cancer cells with advantageous mutations due to genomic instability can be selected and enriched in a Darwinian-type evolutionary process [44]. These mutations are important and recent work has identified point mutations and translocations of the ERα gene (ESR1) in hormone-resistant and metastatic breast cancer [45–47]. Our interest is in cases of reversible resistance that have been observed in experiments and clinical observations. Reversible resistance cannot be explained by genetic mutation and likely occurs on a much shorter timescale. An appealing alternative explanation for this case is the cancer stem cell (CSC) model, which uses the principles of developmental and stem cell biology to explain the multiple cell states and phenotypes as non-genomic sources for Darwinian selection [13]. Increasing evidence of non-genomic heterogeneity has been reported in different cancers, including breast cancer [39]. To determine the efficacy of cancer treatments, there is a significant need to investigate mechanisms of phenotype transitions in cancer cells, and to relate these transitions at the cellular level to the effects of drug selection at the population level.

Considering the three phenotypes of oestrogen (E2) sensitivity in breast cancer cells as an example of non-genomic heterogeneity, we modelled the crosstalk among the three major ERα signalling modes (E2ER, ERM and ERP) and the GFR pathway for cell survival and proliferation. The model, though necessarily high level, provides a qualitative match to the progressive, step-wise phenomenon of shifting E2 sensitivities and changing endocrine responsiveness in breast cancer cells under E2 manipulation.

The complexity of cancer must also be considered. Cancer cells may have as many steady states as there are physiological states, including a spectrum of cell phenotypes from different differentiation stages [48]: CSCs, intermediate progenitor cells and fully differentiated ‘bulk’ cells. Such physiological phenotypes, and the transitions among them, have been widely reported in breast cancer. For example, the epithelial–mesenchymal plasticity of breast cancer creates a broad spectrum of phenotypes, whereby epithelial primary cancers become mesenchymal-like and re-epithelialize into a solid tumour at secondary, metastatic sites [49]. In addition, the recent molecular classification of breast cancers into simple, fixed categories (luminal A, luminal B, HER2+, basal-like, normal-like) is not robust and has neglected the heterogeneous and dynamic properties of cancer [50]. Phenotypic transitions among stem-like, basal and luminal subpopulations within a given breast cancer cell population have been reported in experiments [39]. Unfortunately, current knowledge of the underlying molecular mechanisms for these phenotypic transitions is limited. Nevertheless, the sources of non-genomic heterogeneity in breast cancer cells are now under extensive study and will ultimately enable the creation of enhanced models with more (meta)-stable steady states corresponding to the various physiological states of breast cancer cells.

Mathematical modelling can help understand the emergent behaviours of cancer physiology at different levels [51–56]. For example, several mathematical models have been developed to address the mechanisms by which cancer cells develop resistance to treatment. While most models attribute the source of Darwinian selection to rare mutations [57,58], some have taken non-genetic heterogeneity into account by modelling spontaneous transitions between CSCs and ‘bulk’ cells at the population level [59,60]. In this study, we adopted a multiscale modelling strategy by first creating a biologically plausible, cellular-level model of transitions among key breast cancer phenotypes. We then derived a state transition model for population-level responses to endocrine therapies from the cellular-level model. This strategy enables a systematic understanding of how adaptations to a changing environment can cooperate at both the cellular and population levels to create resistance to oestrogen-withdrawal therapies. Our modelling results indicate that selection pressure can substantially change the shape of the oestrogen response landscape (figure 2*c*,*d*) and change the steady-state distribution of the state-transition model (figure 4*a*,*b*). The transition rate approximation we propose provides a useful bridge between the molecular mechanisms of cellular phenotypic transitions and the phenomenological parameters that have been widely used to describe these transitions in population-based models.

We applied the state transition model to understand the effects of sequential and intermittent treatment that have proved beneficial for breast and other cancers [3,5,7–9,61–63]. Many studies have focused on the benefits of these strategies in reducing side effects. Here we focused on another major advantage, the potential of these strategies to reduce or reverse acquired resistance during treatment. The approach used in this study provides a potential path to investigate treatment strategies and improve treatment efficacy in breast and other cancers. The next step will be to move from a qualitative to a quantitative model in breast cancer cell lines with the goal of making and validating model predictions of the efficacy of therapy protocols. If optimized, time-varying protocols provide significant advantages in reducing resistance and proliferation in cell lines, then the challenging task of moving these results to the clinic can begin.

## Funding statement

This work was supported in part by US National Institutes of Health grant U54-CA149147 (to R.C., J.J.T. and W.T.B.), US National Science Foundation grant DMS-0969417 (to J.X.), and by a fellowship to C.C. provided by the Virginia Polytechnic Institute and State University graduate program in Genetics, Bioinformatics and Computational Biology.

## Acknowledgements

We thank Jin Wang, Lei Zhang, Hang Zhang and Xiaojun Tian for help and/or discussions.

- Received February 26, 2014.
- Accepted April 15, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.