Understanding the origins of resistance to anti-oestrogen drugs is of critical importance to many breast cancer patients. Recent experiments show that knockdown of GRP78, a key gene in the unfolded protein response (UPR), can re-sensitize resistant cells to anti-oestrogens, and overexpression of GRP78 in sensitive cells can cause them to become resistant. These results appear to arise from the operation and interaction of three cellular systems: the UPR, autophagy and apoptosis. To determine whether our current mechanistic understanding of these systems is sufficient to explain the experimental results, we built a mathematical model of the three systems and their interactions. We show that the model is capable of reproducing previously published experimental results and some new data gathered specifically for this paper. The model provides us with a tool to better understand the interactions that bring about anti-oestrogen resistance and the effects of GRP78 on both sensitive and resistant breast cancer cells.
The majority of newly diagnosed breast cancers in women are oestrogen receptor (ERα) positive. These cancers are often treated with anti-oestrogen drugs such as tamoxifen (TAM) or ICI 182,780 (ICI; Faslodex; Fulvestrant), which block the normal functions of ERα or degrade it, resulting in cell cycle arrest or cell death [1,2]. Unfortunately, most patients receiving anti-oestrogen treatment ultimately acquire resistance during treatment . While many molecular components and interactions relevant to anti-oestrogen resistance have been identified experimentally, how these components interact as a system to produce particular phenotypes is still an open question. Moreover, cross talk and feedback among these components makes reasoning about the system difficult and error prone. In this context, mathematical modelling based on experimental data provides a framework for us to analyse oestrogen signalling in a rigorous and consistent manner. The ultimate goal is an integrated model of the oestrogen receptor signalling system that explains the underlying mechanisms of endocrine resistance and enables accurate prediction of the impact of therapies on cell fate. The model may be useful for initial in silico exploration of a wide range of multi-drug therapies to combat and reverse resistance, an exploration that may be impractical to carry out experimentally. Moreover, developing such a model is an effective means for determining gaps and inconsistencies in our current knowledge and suggesting experiments for resolving these problems.
Recent experimental results point to a potentially major role for GRP78, a key player in the unfolded protein response (UPR), in anti-oestrogen resistance . Consistent with previous studies [5,6] that observed increased GRP78 expression in several human breast cancer cell lines, these results showed that knockdown of GRP78 in ICI-resistant cells re-sensitized the cells to ICI, which reduced proliferation due to increased apoptosis. Furthermore, overexpression of GRP78 in an ICI-sensitive cell line conferred resistance to ICI while also upregulating autophagy. These experiments implicate the interaction of three cellular systems, UPR, autophagy and apoptosis, as critical to the role of GRP78 in anti-oestrogen resistance. The experiments also explored key cross-talk mechanisms among these systems by showing that GRP78 modulated mTOR activity (an activator of cell growth and an inhibitor of autophagy) and CHOP expression (a driver of apoptosis). (Full names of proteins, such as mTOR and CHOP, can be found in the electronic supplementary material, table S1.)
Activated by the accumulation of unfolded proteins in the endoplasmic reticulum (EnR) as a result of cellular stresses, the UPR attempts to increase the protein folding and protein degradation capacities of the EnR by upregulating synthesis of chaperones and other proteins [7,8]. A key chaperone in the EnR is GRP78, also known as BiP or HSPA5. Under unstressed conditions, GRP78 binds to EnR stress sensors IRE1, PERK and ATF6, and inhibits their activities . When unfolded proteins accumulate, GRP78 binds to them and attempts to help them fold correctly. GRP78 also participates in the EnR-associated degradation (ERAD) process in which unfolded proteins are retro-translocated across the EnR membrane into the cytoplasm and degraded by proteasomes. The binding of GRP78 to unfolded proteins releases the active forms of IRE1, PERK and ATF6. Activated IRE1, PERK and ATF6 upregulate not only chaperone production but also several effectors of autophagy and apoptosis. Autophagy, a mechanism by which cells degrade unfolded or aggregated proteins and damaged organelles [10–12], reduces stress on the EnR by providing additional energy and raw materials to the cell. If homeostasis cannot be restored in response to the increased stress by means of UPR and autophagy, apoptosis may be triggered. Apoptosis, a form of programmed cell death, is an orderly process of cell death that can be initiated by either internal or external signals.
Several recent studies report on the complex interplay between UPR, autophagy and apoptosis that decides cell fate in response to various drug treatments [4,13–23]. Anti-oestrogens or other drug treatments can activate either cyto-protective or cyto-destructive UPR, and autophagy has been shown to aid in the cyto-protective role of the UPR [21–24], whereas cyto-destructive UPR leads to apoptosis [18,25].
The aim of this study is to determine whether our current mechanistic understanding of relevant cellular interactions is sufficient to explain the experimental results concerning the role of GRP78 in anti-oestrogen resistance reported in Cook et al. . Our approach to modelling the interaction of these three cellular systems in response to GRP78 manipulation and ICI treatment starts with building models of the individual subsystems using data from the literature to arrive at baseline models that are not based on breast cancer cells. We then add cross-talk pathways to form an integrated model of the signalling network in breast cancer cells , and tune the model based on experimental data from two cell lines that reflect the clinical development of acquired endocrine resistance: LCC1, an oestrogen-independent but anti-oestrogen-sensitive line derived from the standard ERα-positive breast cancer cell line MCF-7 ; and LCC9, an oestrogen-independent and ICI-resistant line derived from LCC1 . The mathematical model replicates the experimentally observed behaviour of the LCC1 and LCC9 cell lines when subjected to GRP78 knockdown or overexpression alone, and when combined with ICI treatment.
2. Model development
As shown in figure 1, the signalling network integrates three cellular systems: the UPR, apoptosis and autophagy. Each of these systems is modelled at a different level of detail, reflecting both differences in the degree of molecular knowledge about each system and the needs of our integrated model. For example, the molecular mechanism associated with the UPR is known in much greater detail than the mechanism of autophagy, and the molecular details of UPR interactions are necessary for understanding the influence of GRP78. By contrast, while the molecular details of the intrinsic apoptosis pathway are fairly well known, only the initial switch-like triggering of apoptosis is needed for the experimental results we are considering.
Similarly, while there are many known pathways that produce interactions among the cellular systems in our model, we have only included those paths that appear to be the dominant ones for the experiments we consider here. In the future, as the model is used to explain additional experimental data, currently excluded pathways may become important and need to be added to the model. A complete list of the equations and parameters for our model, along with simulation code, is included in the electronic supplementary material.
2.1. Unfolded protein response module
Our representation of the UPR module is adapted from a model by Verdugo and co-workers . The model captures the binding of chaperones to unfolded proteins and the resulting reduction in unfolded proteins through either correct folding or degradation (ERAD; figure 1). The binding of chaperones to unfolded proteins releases the three transmembrane sensors IRE1, PERK and ATF6 from the repressing effects of chaperones  and activates the three arms of the UPR. (i) Activated IRE1 converts unspliced mRNA (mXBP1u) for X-box-binding protein 1 to a spliced form (mXBP1s), which produces a transcription factor (XBP1s protein) that upregulates GRP78. (ii) Activated PERK phosphorylates eukaryotic initiation factor-2 (eIF2), which inhibits the translation of most proteins, thereby reducing the load of unfolded proteins entering the EnR. However, the translation of a few proteins, including ATF4 and CHOP, is increased by the phosphorylation of eIF2. ATF4 is a transcription factor that can increase CHOP expression, and CHOP regulates the transcription of proteins responsible for apoptosis . (iii) Activated ATF6 is transported to the Golgi apparatus, where it is cleaved to produce a transcription factor that upregulates synthesis of XBP1 mRNA. ATF4, XBP1s and cleaved ATF6 also upregulate numerous EnR chaperones, which increase the protein-folding capacity of the EnR.
We assume that the association of UPR sensors with GRP78 is not the sole regulatory mechanism of IRE1, PERK and ATF6 activity, since these sensors can still upregulate the EnR stress response in cells manipulated by either GRP78 knockdown or disrupted binding of GRP78 with the sensors [4,29]. These results suggest that, in the absence of GRP78, alternate inhibitory mechanisms (e.g. other chaperones) release the sensors for further activation upon additional EnR stress (reviewed in Ma & Hendershot ). While numerous chaperones are likely to play a role in sensor inhibition, such as BI-1 [31,32], AIP1  and RACK1  for IRE1, and WFS1  for ATF6, for simplicity of the model, we lump the effects of these other chaperones into the effect of a single chaperone, GRP94. While not redundant to GRP78, GRP94 plays an important role in EnR stress control [36–38]. Furthermore, GRP78 knockdown strongly induces GRP94 expression [4,39], indicating that cells can use an alternative mechanism to increase protein-folding capacity. To account for the effect of other chaperones, we modify Verdugo's UPR model by including GRP94 as an alternative stress regulator. While other chaperones play an essential role in EnR homeostasis, their modes of action on UPR sensors are still obscure relative to the well-characterized role of GRP78. Hence, we assume that GRP94 controls the UPR arms in the same way as GRP78, and modify the equations accordingly. We assume a constant level of GRP94 protein, ignoring its upregulation by EnR stress, since its enhanced expression under GRP78 knockdown does not relieve the stress, as discussed in Cook et al. , and so its upregulation may not be significant.
2.2. Autophagy module
Autophagy is principally an adaptive response to protect cells from stress, but, in certain settings, it can lead to cell death [40,41]. The process is controlled by a set of evolutionarily conserved autophagy-related proteins (ATG proteins). ATG6/Beclin-1 (BECN1), ATG8 (also known as microtubule-associated protein 1 light chain 3, LC3), ATG5, ATG12, ATG9 and ATG13 are essential for autophagosome formation in mammalian cells [42–44]. We model only the initial step of autophagy, namely the initiation of an isolation membrane, which is subsequently elongated to form the autophagosome, a vesicle enclosing the cellular material to be recycled. Autophagosomes then fuse with lysosomes, and their contents are degraded. BECN1 and ATG13 are indispensable for the formation of the isolation membrane. Under unstressed conditions, BCL2EnR (BCL2 in the EnR) binds to BECN1 and inhibits its activity [20,44,45]. Phosphorylation of BCL2EnR by activated JNK releases BECN1 from BCL2EnR, which then leads to an induction of autophagy. The active mammalian target of rapamycin complex (mTORC1*) negatively regulates autophagy by suppressing the activity of the ATG1 : ATG13 : ATG17 complex [44,46,47]. mTORC1 is regulated by a variety of upstream components such as AMPK, TSC1/2, AKT and intracellular Ca2+, depending upon the type of stress [4,15,48–50]. The transport of Ca2+ from the EnR to the cytoplasm, which plays a key role in the model, is regulated by the inositol-1,4,5-triphosphate receptor (IP3R), a Ca2+ channel in the EnR membrane .
Our model of autophagy initiation, based on a model by Tavassoly and co-workers , comprises the following steps (figure 1): (i) phosphorylation of BCL2EnR by phosphorylated JNK, and subsequent release of BECN1 and IP3R from inhibition by BCL2EnR, (ii) deactivation of mTORC1 by cytoplasmic Ca2+ and GRP78 , relieving its repressive effect on ATG13 and inducing expression of ATG9, and (iii) formation of autophagosomes by BECN1, ATG13 and ATG9 interaction. We have added ATG9 to Tavassoly's model because experiments show increased expression of ATG9 and increased autophagy upon GRP78 overexpression. Without ATG9, overexpression of GRP78 in the model leads to repression of IRE1 activity, inhibition of JNK phosphorylation, decreased binding of BECN1 to ATG13 and ultimately to downregulation of autophagy. Including ATG9 provides an alternative pathway that compensates for this downregulation.
2.3. Apoptosis module
We consider the intrinsic apoptosis pathway in which mitochondrial outer membrane permeabilization (MOMP) represents the key decision process and point of no return . Ignoring downstream details, we model only the switch-like MOMP event by considering three classes of proteins: BCL2mito (anti-apoptotic BCL2 family members in mitochondria), BH3 (BH3-only BCL2 family members) and BAX (pro-apoptotic BCL2 family members directly involved in MOMP). Under non-apoptotic conditions, BCL2mito binds and sequesters BAX in an inactive form. The intrinsic apoptosis pathway responds to elevated intracellular stresses by activating BH3 proteins, which release BAX from the inhibitory effect of BCL2mito. Released BAX self-oligomerizes and forms pores in the mitochondrial outer membrane, releasing the mitochondrial intermembrane protein cytochrome c (Cyt-C) into the cytosol [52,53]. Cyt-C activates pro-caspases that subsequently activate effector caspases . The model includes interactions between BCL2mito, BH3 and BAX as shown in figure 1 [26,54]. In the model, initiation of apoptosis is signalled by the rapid increase in free BAX, whereas, in the experiments, it was quantified in terms of the levels of both cleaved CASPASE-7 (a downstream effector caspase) and annexin V staining.
2.4. Integration of unfolded protein response, autophagy and apoptosis modules
EnR stress is a potent trigger of autophagy [16,55]. Autophagy is usually considered a pro-survival response, protecting cells against starvation, infectious agents and treatment with anti-cancer drugs. However, under certain circumstances, autophagy can promote death [40,41]. When EnR stress is extensive or sustained, and the UPR fails to restore normal functioning, apoptosis is triggered . Several studies have shown that apoptosis and autophagy may be interconnected in some settings [56–60].
IRE1 and increased cytoplasmic Ca2+ have been implicated as inducers of autophagy linked to EnR stress in mammalian cells [15,16]. Activation of the IRE1 branch of the UPR can lead to phosphorylation and activation of JNK [16,61,62], and activated JNK phosphorylating BCL2EnR leads to the dissociation of BCL2EnR from BECN1 and activation of autophagy . EnR stress can also lead to release of Ca2+ from the EnR into the cytosol through the IP3R channel, which in turn can activate autophagy by inhibiting mTORC1 activity [15,63]. BCL2EnR has been proposed as a regulator of IP3R function, and thus the cytoplasmic Ca2+ level [15,45]. In particular, phosphorylation of BCL2EnR by phosphorylated JNK renders the IP3R active and increases the cytoplasmic Ca2+ concentration . These pathways are shown in figure 1.
In addition, Cook et al.  proposed that GRP78 inhibits mTOR phosphorylation and decreases mTORC1 activity by activating AMPK and phosphorylating TSC2 (also reviewed in Cook & Clarke ). For simplicity, the details of mTOR phosphorylation and mTORC1 formation are modelled as a simple reaction that activates or deactivates the mTORC1 complex. Thus, the effect of GRP78 on mTORC1 has been modelled as a direct effect, excluding the unknown intermediate reactions between GRP78, AMPK and TSC2. To account for the pro-survival effect of autophagy, we add an autophagy-dependent negative feedback into the equation for cellular stress.
The UPR is usually linked to apoptosis through its PERK signalling branch, which downregulates the expression of anti-apoptotic BCL2mito and upregulates pro-apoptotic BH3 proteins [66,67]. In the model, this regulation is accomplished via increased expression of CHOP, which has been shown to promote apoptosis in several contexts [18,25,66,67]. The UPR can also signal to apoptosis via an IRE1–JNK pathway. We have not explicitly included this pathway in our initial modelling, lumping all UPR signalling to apoptosis into the CHOP pathway.
In addition, increased cytoplasmic Ca2+ can stimulate the production of reactive oxygen species (ROS) in mitochondria, impair electron transport and oxidative phosphorylation and increase expression of pro-apoptotic BH3 proteins. In the model, we lump these effects of cytoplasmic Ca2+ on apoptosis into its action on BH3 protein expression. This role for cytoplasmic Ca2+ in triggering apoptosis via the UPR is consistent with a study by Ouyang et al. , who observed that increasing GRP78 protects astrocytes from ischaemic injury by decreasing the net flux of Ca2+ from the EnR to the mitochondria.
2.5. ICI treatment
ICI treatment increases ERα aggregation in the cytoplasm, which generates a weak signal that activates the UPR . Furthermore, Cook et al.  demonstrated that ICI treatment induces a low level of autophagy in breast cancer cells, consistent with its modest effect on the UPR. Therefore, in the model, we consider a weak activation of UPR by ICI that comes through a parameter ‘stress’, which represents the rate of production of unfolded proteins.
In response to ICI treatment of anti-oestrogen-sensitive breast cancer cells, the major apoptotic stimulus is through ERα degradation not UPR activation. Expression of the anti-apoptotic protein BCL2mito is promoted by ERα , and since ICI degrades ERα, ICI therapy results in a decreased level of BCL2mito, favouring apoptosis . In the model, the BCL2mito level is controlled by ERα and NFκB. We do not model details of ERα signalling but merely allow for its synthesis and degradation, including an ICI-dependent degradation term.
The addition of ICI to anti-oestrogen-resistant cells that have had GRP78 knocked down increases mTORC1 activity and decreases autophagy . The exact mechanism of mTORC1 activation is unknown. However, it is known that ICI can induce AKT or ERK1/2 phosphorylation in ERα-negative breast cancer cells  or in ERα-positive and HER2-abundant breast cancer cells , which can result in mTOR phosphorylation and enhanced mTORC1 activity. This ICI-mediated activation of AKT/ERK depends on ERα36, a membrane-bound ERα  that is not degraded by ICI. However, given the complex bidirectional cross talk between EnR stress and mTOR signalling , we cannot rule out other mechanisms involving downstream effectors of UPR. In the model, we add a phenomenological term to account for the effect of ICI on mTORC1.
The deterministic model we have developed captures the average behaviour of one cell in a population of N cells. To model the expansion or contraction of the population, we use a simple birth–death equation of the form 2.1where the specific birth and death rates of the population, kB and kD, are functions of the signalling network that determines cellular responses to stress. In particular, phosphorylated (active) mTOR, through the active protein complex mTORC1*, regulates global rates of protein synthesis. Since the rate of cell division must match the rate of cell growth, we write the specific birth rate kB as a function of activated mTORC1 in the form 2.2where T0 is the intrinsic doubling time of the cell culture (28 h). With our choice of parameter values (w1 = 1.9, wmTORC1 = 0.07), kB increases smoothly from 0.0215 h−1 at [mTORC1*] = 0 to 0.0247 h−1 at [mTORC1*] (see the electronic supplementary material, figure S1(a)).
The specific death rate in experiments and in the model is determined by the rate of apoptosis. In our model, the rate of apoptosis is determined by the relative levels of BCL2mito and BH3 proteins. Individual cells in a population will have varying levels of these proteins owing to stochastic effects within the cells. Cells with higher levels of BCL2mito or lower levels of BH3 will be more immune to apoptotic stimuli, whereas cells with less BCL2mito or more BH3 will be more likely to die in response to such stimuli. For this reason, we write the specific death rate kD as a function of BCL2mito and BH3 in the form 2.3where is the specific death rate of apoptotic cells (when [BH3] [BCL2mito]). With our choice of parameter values (w2 = 0.1, wBCL2 = −0.07, wBH3 = 0.05), when [BH3] = 1.4 × [BCL2mito] + 12. For smaller values of [BH3], kD drops rapidly to approximately 10−3 h−1 (see the electronic supplementary material, figure S1(b)).
To compute the expansion of a cell population over a specified time period, 0 ≤ t ≤ Tend, we first compute the steady-state values of active mTORC1, BCL2 and BH3 from the model under specific experimental conditions. Then, we compute kB and kD according to equations (2.2) and (2.3). Finally, equation (2.1) implies that N(Tend) = N(0) · exp[(kB − kD) · Tend].
2.7. Matching experimental data
Our baseline model refers to LCC1 cells. To model LCC9 cells, we make a small number of parameter changes to account for the biological differences observed between these two cell lines. In anti-oestrogen-resistant LCC9 cells, growth factor signalling replaces ERα signalling as the driver of proliferation [75–77]. Also, compared with LCC1 cells, LCC9 cells have higher NFκB activity  and lower levels of expression of ERα . Thus, for modelling LCC9 cells, the synthesis rate of the ERα is decreased and the amount of active NFκB is increased relative to the LCC1 model. In addition, to account for increased GRP78 levels in LCC9 cells, we assume a higher basal level of stress and hence of unfolded protein production than in LCC1 cells. This assumption is consistent with the fact that LCC9 cells have higher basal autophagy than LCC1 cells .
In the experimental protocol, transient transfection was performed 24 h before ICI treatment, and ICI treatment was continued for 6 days before final measurements were made . To reproduce this protocol in our simulations, the initial conditions (at t = −24 h) for the simulation were chosen to correspond to a steady state of the model, knockdowns or overexpressions were initiated at −24 h and ICI was added at t = 0. The model was then simulated for an additional 144 h (6 days) to allow direct comparison between the model endpoint and the experimental measurements.
The subsystem models and cross-talk mechanisms discussed above were merged into an overall model and the parameters were adjusted to reproduce the experimental measurements at 144 h reported by Cook et al. . In addition, some time-course data were gathered specifically for this study, to ensure that the model captures key transient characteristics of the system in addition to capturing its end-to-end performance. We compare model simulations with observations on four cell types: LCC1, LCC1/GRP78+, LCC9 and LCC9/GRP78−. After showing that the model is consistent with observed cell proliferation for all cases, we consider each case individually and show that the model captures key internal states of the system related to apoptosis, growth (mTORC1* level) and autophagy.
When LCC1 cells are treated with increasing doses of ICI, they exhibit a marked decrease in relative proliferation at 6 days (figure 2a). By contrast, LCC9 cells show no change in proliferation as a function of ICI dose (figure 2b). Figure 2 also shows that the sensitivity of LCC1 cells and the resistance of LCC9 cells can be reversed by GRP78 overexpression and GRP78 knockdown, respectively. The end-to-end performance of our model for these four cases shows good agreement with the experimental results in figure 2.
For the case of LCC1 cells treated with ICI, the reduction in proliferation appears mainly to be the result of increased levels of apoptosis, as illustrated by the marked increase in the number of cells positive for annexin V staining in response to 100 nM ICI (figure 3a). A simulation of BAXmT (the active form of BAX in our model) in response to treatment with 100 nM ICI (figure 3b) indicates that an ‘average’ cell would commit to apoptosis near the 24 h time point. The apoptotic response is due mainly to the degradation of ERα by ICI and the resulting decrease of anti-apoptotic BCL2mito. In our deterministic model, we simulate a single cell, and for specified initial conditions, the cell either commits to apoptosis or it does not. In a population of cells, however, a specified treatment will cause some cells to commit to apoptosis and others to live, owing (presumably) to the fact that protein levels in individual cells are distributed over a range of values [80,81]. Thus, cells with lower levels of BCL2mito are more likely to commit apoptosis, whereas those with higher levels are less likely to do so. The model captures this effect in the formula for the specific death rate, kD, in the proliferation equation.
In terms of deterministic simulations, movement of apoptosis initiation earlier in time indicates that even more BCL2mito would be required to avoid apoptosis for this case; thus, we would expect a larger fraction of cells in a population to die than in a simulation in which apoptosis was initiated later. Figure 3b shows that the initiation of apoptosis in the model does indeed move earlier in response to increasing ICI dose. Although ICI decreases the proliferation rate mainly by inducing apoptosis, the experimentally observed decrease in phospho-mTOR after ICI treatment (figure 4a) may also contribute to the reduced proliferation. The model captures the effect of ICI on mTOR through the reduction of active mTORC1 (figure 4b), and the effect on growth through the dependence of the specific birth rate, kB, on mTORC1*. We assume that mTORC1 activity, which we model directly, reflects the active (phosphorylated) mTOR level. For the case of treatment with 100 nM ICI, our model predicts a decrease in proliferation of 42 per cent (figure 2a), with 37 percentage points of the decrease coming from increased apoptosis and 5 percentage points from decreased growth, showing that apoptosis is the dominant mechanism.
ICI treatment modestly increases autophagy in LCC1 cells in experiments (figure 5a) and model simulations (figure 5b). Experiments blocking autophagy via ATG5 knockdown enhanced the killing effect of ICI , showing that the ICI-induced autophagic response is pro-survival. We simulate these experiments by knocking down ATG9 in our model, which causes apoptosis to initiate earlier in time (figure 3b), indicating that a larger fraction of cells will die. The reduction in autophagy upon ATG9 knockdown decreases the negative feedback to stress and thereby increases the ICI-induced stress level. The increased stress intensifies UPR activation resulting in increased CHOP expression, which drives apoptosis earlier.
For the case of LCC1 cells overexpressing GRP78, ICI treatment does not induce apoptosis, as shown in figure 3a. To produce the corresponding simulation (figure 3b), GRP78 was overexpressed by 3.5-fold, in accord with the experimentally observed level in fig. 2d of Cook et al. . We see no triggering of apoptosis and thus expect no substantial decrease in proliferation due to cell death. While ICI treatment decreases ERα levels and thereby BCL2mito, apoptosis in this simulation is averted by GRP78 overexpression through three mechanisms: (i) a decrease in pro-apoptotic signalling from PERK to CHOP, (ii) an increase in autophagy via GRP78-mediated mTORC1 inactivation and consequent reduction in UPR activity, and (iii) a decrease in pro-apoptotic signalling from IRE1 to cytoplasmic Ca2+ via phospho-JNK. To determine the relative strengths of these mechanisms in the model, we can hold the levels of CHOP, autophagy or Ca2+ equal to their levels in ICI-treated LCC1 cells at 144 h and examine how this changes proliferation for the LCC1/GRP78+ case. For treatment with 100 nM ICI, the proliferation of LCC1/GRP78+ cells decreases by 51, 28 or 25 per cent when CHOP, autophagy or Ca2+, respectively, are held at LCC1 levels, indicating that a decrease in pro-apoptotic signalling from PERK to CHOP is the dominant mechanism.
Figure 4a shows that GRP78 overexpression dramatically suppresses mTOR activity in experiments. ICI treatment, either alone or in combination with GRP78 overexpression, has no additional effect. Model simulations (figure 4b) exhibit the same qualitative response as experiments. In the model, the reduction in proliferation owing to GRP78 overexpression (observed in figure 2a) is due to a slowdown in cell growth captured by the dependence of the birth rate on active mTORC1. Biologically, GRP78-induced dephosphorylation of mTOR, and hence decreased mTORC1 activity, causes the phosphorylation and inactivation of initiation factors eIF4E and eIF2E, which reduces the rate of protein synthesis in cells .
Because mTORC1* significantly inhibits autophagy, both experimental and simulation results show that autophagy is markedly upregulated by GRP78 overexpression (figure 5a,b). Experiments further show that ATG5 knockdown restores sensitivity of LCC1/GRP78+ cells to ICI, suggesting a significant role for autophagy in the resistance of LCC1/GRP78+ (results not shown). Similarly, BAXmT simulations (figure 3b) show that overexpression of GRP78 eliminates apoptosis, whereas knocking down autophagy through ATG9 restores the initiation of apoptosis in response to ICI treatment.
Experimental fold-changes at the 6 day time point due to GRP78 overexpression were reported for mTORC1* (down 6.9-fold), CHOP (down 3.0-fold), ATG9 (up 3.0-fold) and total BCL2 (up 2.0-fold) in Cook et al. . Model simulation results are in good agreement: mTORC1* (down 7.0-fold), CHOP (down 3.1-fold), ATG9 (up 2.6-fold) and total BCL2mito (up 2.1-fold). In the model, total BCL2mito was taken as the sum of BCL2mito, BCL2 : BH3 and BAXm : BCL2 levels.
Normal LCC9 cells are insensitive to ICI treatment (figure 2b), as seen by the lack of apoptosis as measured by annexin V staining (figure 3c). Simulation of an LCC9 cell with ICI treatment shows no activation of BAX, and hence no apoptosis, over 6 days (figure 3d). LCC9 cells behave this way in our simulation because their BCL2 level is elevated by their increased NFκB activity. (NFκB directly binds the promoters and induces the expression of several anti-apoptotic proteins, including BCL-2 family members .) LCC9 cells are also less sensitive to degradation of ERα by ICI since they have less ERα to begin with. Figure 6d shows that the mTORC1* level in our model of LCC9 cells does not change significantly in response to ICI, so cell growth rate is not expected to change significantly. The combination of no apoptosis and unchanged growth rate renders our simulated LCC9 cells insensitive to ICI. ICI treatment produces a small increase in autophagy in both experiments and simulations (figure 5c,d). As might be expected, knocking down autophagy via BECN1 has little effect on the insensitivity of LCC9 cells to ICI in the experiments or simulations (figure 2b), although high doses of ICI do reduce proliferation somewhat.
When GRP78 is knocked down in LCC9 cells, proliferation decreases significantly (approx. 25%) and treatment with ICI causes an intensifying decrease with increasing dose (figure 2b). These effects are due to increasing levels of apoptosis, as indicated by the experimentally observed increase in annexin V staining (figure 3c). Model simulations of these experiments (figure 3d) show that BAX is activated when GRP78 is knocked down and adding ICI moves the activation earlier in time; thus, a larger fraction of cells in a population will die. The increases in apoptosis come from four sources: (i) knocking down GRP78 in the model increases activation of the UPR, which in turn provides pro-apoptotic signalling via increased CHOP expression and increased cytoplasmic Ca2+ levels; (ii) ICI dramatically increases mTORC1 activation (figure 6d), inhibiting autophagy (figure 5d) and causing increased UPR activation and pro-apoptotic signalling; (iii) ICI reduces the level of ERα, via enhanced degradation, leading to a decrease in the level of BCL2mito that increases the likelihood of apoptosis; and (iv) ICI increases EnR stress and causes a further activation of apoptotic signalling.
The experimental results in figure 6d indicate no dramatic changes in mTORC1 activity when GRP78 is knocked down, so we expect that the changes in proliferation observed in figure 2b are due to changing levels of apoptosis. The experimental results also show that adding ICI to LCC9/GRP78(−) cells significantly intensifies the effect of GRP78 knockdown on mTORC1 activation. Simulations of the model in figure 6d are in good agreement with these results. Simulations of normal LCC9 cells with ICI treatment (figure 6d) show little change in mTORC1 activity because the initial increase in mTORC1 activity by ICI is offset by ICI activation of UPR, resulting in an increased GRP78 expression and cytosolic Ca2+ level that both inhibit mTORC1 activity. The dramatic increase when ICI is combined with GRP78 knockdown is a result of the direct action of ICI on mTORC1 activation in the model combined with the saturation of the mTORC1 deactivation reaction.
Experimentally, we also find that GRP78 knockdown has little effect on autophagy after 6 days (figure 5c), whereas the model shows a small decrease from the baseline value (figure 5d). The addition of ICI in the presence of a GRP78 knockdown in LCC9 cells significantly reduces the extent of autophagy in both experiments and simulations (figure 5c,d). This reduction in the model is a consequence of ICI's role in promoting mTORC1 activation. Furthermore, knocking down both autophagy (via BECN1) and GRP78 decreases proliferation in experiments and simulations (figure 2b).
We have collected time-course data for four key proteins, GRP78, phospho-JNK, CHOP and mTORC1*, for LCC9 cells treated with GRP78 siRNA alone and in combination with ICI (symbols, figure 6). The data show the fold change of each protein, compared with its unstressed level, at several time points. Our simulations capture the experimental data quite well (lines, figure 6). The experimental data in figure 6a show that GRP78 siRNA successfully decreases GRP78 protein expression by 80 per cent within 72 h. This decrease activates the UPR and ultimately causes JNK phosphorylation and increases CHOP expression (figure 6b,c). Owing to the inverse relationship between GRP78 expression and mTORC1 activity, GRP78 knockdown stimulates mTORC1 activation (figure 6d). ICI treatment of GRP78(−) cells does not have a major effect on GRP78 expression (figure 6a), but it increases phospho-JNK, CHOP expression and mTORC1 activity significantly compared with GRP78 knockdown alone (figure 6b–d).
Our model provides a clear, consistent and quantitatively accurate account of how we think apoptosis, autophagy and UPR interact to produce both sensitivity and resistance to anti-oestrogen therapy under various conditions. The model allows us to characterize the relative importance of the multiple pathways that create the observed phenotypes and demonstrates that the mechanisms we have chosen to model are sufficient to explain the experimental results under consideration. However, this success does not prove that other pathways are unimportant to the current experimental results or will not need to be considered as we move forward with modelling other cellular responses to anti-oestrogen therapy.
An interesting aspect of the model, and the experiments, is that, whether GRP78 is overexpressed (in LCC1 cells) or knocked down (in LCC9 cells), proliferation in both cases decreases when there is no ICI treatment. This is due to the nonlinear nature of the model. If GRP78 is sufficiently overexpressed, proliferation ultimately decreases owing to inactivation of mTOR. Similarly, if GRP78 is sufficiently knocked down, proliferation ultimately decreases owing to increasing apoptosis. Somewhere between these two extremes, there is an optimal level of GRP78 at which proliferation is a maximum. For cells not at this maximum level, if a small decrease in GRP78 causes proliferation to increase, then a small increase in GRP78 will cause proliferation to decrease. This is the case for LCC9 cells. But for large changes in GRP78, proliferation ultimately decreases regardless of whether GRP78 is increased or decreased.
In the model, GRP78 overexpression creates resistance by suppressing pro-apoptotic signalling, using the three mechanisms pointed out in §3. However, contributions to resistance from other mechanisms under the direct control of GRP78 cannot be excluded. For example, GRP78 can suppress EnR stress-induced apoptosis by sequestering caspase-7 and inhibiting its activation [84,85]. GRP78 can also inhibit pro-apoptotic protein BIK . There may also be a contribution from growth factor receptor (GFR) signalling to GRP78 overexpression-induced resistance, since membrane-bound GRP78 can activate RAS–MAPK and PI3K–AKT pathways in other cells [87,88]. The effect of growth factor signalling may not be significant in our cells, as Cook et al.  observed no change in anti-oestrogen resistance when cell surface signalling was blocked in LCC9 cells. Measurements of AKT and ERK activity in LCC1/GRP78(+) cells may clarify whether GFR signalling is important or not.
To model LCC9 cells, we assumed a higher basal rate of unfolded protein production than in LCC1 cells. Persistent ROS stress or persistent production of unfolded proteins may induce adaptive stress responses, enabling cells to maintain viability . Since LCC9 cells were obtained from LCC1 cells by a stepwise selection using ICI, we speculate that LCC9 cells are adapted to survive higher stress levels than their parental LCC1 cells. Furthermore, LCC9 cells have a higher basal level of autophagy than LCC1 cells, indicating higher basal stress. The increased GRP78 mRNA and protein levels in LCC9 cells may be a consequence of mechanisms other than just increased basal stress, however. Dai et al.  found that PI3K/AKT promotes increased accumulation of GRP78 in HEK293 cells by increasing protein stability and not by increasing the GRP78 mRNA level. This may be the case in LCC9 cells as well, since they have high PI3K/AKT activity . Other UPR-independent mechanisms that can result in increased GRP78 mRNA and protein levels include IGF-1R/PI3K/AKT/mTORC1 signalling  and Raf1–MEK–ERK signalling , observed in mouse embryonic fibroblast cells and in human gastric tumour cells, respectively. It would be interesting to know whether the increased level of GRP78 is a consequence of epigenetic changes or is actively maintained to deal with higher rates of unfolded proteins.
In the model, GRP78 knockdown restores anti-oestrogen sensitivity in LCC9 cells by increased apoptotic signalling through CHOP and Ca2+. However, other mechanisms may be involved. As mentioned above, GRP78 can sequester caspases and pro-apoptotic proteins, so GRP78 knockdown may release these proteins and enhance apoptotic signalling. In addition, EnR stress-induced signalling to apoptosis via CASPASE-4/12 [94,95] or the IRE1–JNK pathway may be important. GRP78 knockdown is also known to induce apoptosis by blocking AKT activation in colon, PTEN-null leukaemia and prostate cancer cell lines [88,96,97]. Hence, it will be important to measure the effect of knocking down GRP78 on AKT signalling in LCC9 cells.
By providing a temporal interaction mechanism capable of explaining key aspects of anti-oestrogen resistance, this model is a step towards a better understanding of an important clinical problem affecting a large number of breast cancer patients. Accurate prediction, however, will require a model that is based on a much larger set of experimental results. Future work will improve the model by considering data involving knockdowns and overexpression of other key resistance players, such as XBP1, NFκB and ERα. In addition, we plan to extend the model to include relevant pathways from growth factor signalling.
5. Experimental material and methods
The following materials were obtained as indicated: ICI 182,780 (Tocris Bioscience, Ellisville, MO); improved minimal essential medium (IMEM; Gibco Invitrogen BRL, Carlsbad, CA); bovine calf charcoal stripped serum (Equitech-Bio Inc., Kerrville, TX); Lipofectamine RNAiMax reagent (Invitrogen); GRP78 siRNA (On-Target plus SMART pool (consisting of three different siRNA for the same target), ThermoScientific Dharmacon, Lafayette, CO); antibodies were obtained from the following sources: GRP78, CHOP, JNK, phospho-JNK and mTORC1 (Cell Signalling Technology); β-actin and polyclonal and horseradish peroxidase (HRP)-conjugated secondary antibodies (Santa Cruz Biotechnology, Santa Cruz, CA).
5.2. Cell culture
MCF7/LCC9 (LCC9) breast carcinoma cells were grown in phenol-red-free IMEM medium containing 5 per cent charcoal-treated calf serum. Cells were grown at 37°C under a humidified 5 per cent CO2/95 per cent air atmosphere. Cells were plated in a six-well plate and transfected with control or GRP78 siRNA. The next day, the medium was replaced and cells were treated with vehicle or 100 nM ICI for 0, 24 or 48 h before protein was harvested.
5.3. Western blot hybridization
Treated cell monolayers were solubilized in radioimmunoprecipitation assay lysis buffer. Protein concentrations of cell lysates were measured using a standard bicinchoninic acid assay. Proteins were size fractionated by polyacrylamide gel electrophoresis and transferred to nitrocellulose membrane. Non-specific binding was blocked by incubation with Blotto (Tris-buffered saline containing 5% powdered milk and 1% Triton X-100). Membranes were incubated overnight at 4°C with primary antibodies. The next day, membranes were washed and then incubated with polyclonal HRP-conjugated secondary antibodies (1 : 2000) for 1 h at room temperature. Immunoreactive products were visualized and quantified by densitometry using the ImageJ digital densitometry software (http://rsbweb.nih.gov/ij/). Protein loading control was determined by incubation of stripped membranes with a monoclonal antibody to β-actin (1 : 1000).
One contribution of 11 to a Theme Issue ‘Integrated cancer biology models’.
- © 2013 The Author(s) Published by the Royal Society. All rights reserved.