Stochastic multi-scale models of competition within heterogeneous cellular populations: Simulation methods and mean-field analysis

We propose a modelling framework to analyse the stochastic behaviour of heterogeneous, multi-scale cellular populations. We illustrate our methodology with a particular example in which we study a population with an oxygen-regulated proliferation rate. Our formulation is based on an age-dependent stochastic process. Cells within the population are characterised by their age (i.e. time elapsed since they were born). The age-dependent (oxygen-regulated) birth rate is given by a stochastic model of oxygen-dependent cell cycle progression. Once the birth rate is determined, we formulate an age-dependent birth-and-death process, which dictates the time evolution of the cell population. The population is under a feedback loop which controls its steady state size (carrying capacity): cells consume oxygen which in turn fuels cell proliferation. We show that our stochastic model of cell cycle progression allows for heterogeneity within the cell population induced by stochastic effects. Such heterogeneous behaviour is reflected in variations in the proliferation rate. Within this set-up, we have established three main results. First, we have shown that the age to the G1/S transition, which essentially determines the birth rate, exhibits a remarkably simple scaling behaviour. Besides the fact that this simple behaviour emerges from a rather complex model, this allows for a huge simplification of our numerical methodology. A further result is the observation that heterogeneous populations undergo an internal process of quasi-neutral competition. Finally, we investigated the effects of cell-cycle-phase dependent therapies (such as radiation therapy) on heterogeneous populations. In particular, we have studied the case in which the population contains a quiescent sub-population. Our mean-field analysis and numerical simulations confirm that, if the survival fraction of the therapy is too high, rescue of the quiescent population occurs. This gives rise to emergence of resistance to therapy since the rescued population is less sensitive to therapy.

for phenotypes to increase their robustness as time progresses). These properties are exploited by tumours to increase their proliferative potential and resist to therapies (Kitano, 2004). In addition to complex, non-linear interactions within individual cells, there exist intricate interactions between different components of the biological systems at all levels: from complex signalling pathways and gene regulatory networks to complex non-local effects where perturbations at whole-tissue level induce changes at the level of the intra-cellular pathways of individual cells (Alarcón et al., 2005;Ribba et al., 2006;Macklin et al., 2009;Osborne et al., 2010;Deisboeck et al., 2011;Powathil et al., 2013;Jagiella et al., 2016). These and other factors contribute towards a highly complex dynamics in biological tissues which is an emergent property of all the layers of complexity involved.
In parallel to the model development, algorithms and analytic methods are being developed to allow for more efficient analysis and simulation of such models (Alarcón, 2014;Spill et al., 2015;de la Cruz et al., 2015;. In the case of cancer biology, the multi-scale interactions of intracellular changes at the genetic or molecular pathway level and tissue-level heterogeneity can conspire to generate unfortunate effects such as resistance to therapy (Merlo et al., 2006;Gillies et al., 2012;Greaves and Maley, 2012;Chisholm et al., 2015;Asatryan and Komarova, 2016). Heterogeneity plays a major role in the emergence of drug resistance within tumours and can be of diverse types. There is heterogeneity in cell types due to increased gene mutation rate as a consequence of genomic instability and other factors (Merlo et al., 2006;Greaves and Maley, 2012;Chisholm et al., 2015;Asatryan and Komarova, 2016). Heterogeneity can also be caused by the complexity of the tumour microenvironment (Alarcón et al., 2003;Gillies et al., 2012;Chisholm et al., 2015), in which diverse factors such as tumoural or immune cells (Kalluri andZeisberg, 2006, Grivennikov et al., 2010), or the extracellular matrix and its physical properties , strongly influence cancer cell behaviour. Note that hypoxia is also known to change the tumour microenvironment . In either case, heterogeneity within the tumour creates the necessary conditions for resistant varieties to emerge and be selected upon the administration of a given therapy.
The main aim of this paper is to analyse the properties of heterogeneous populations under the effects of fluctuations both within the intracellular pathways which regulate (individual) cell behaviour and those associated to intrinsic randomness due to finite size of the population. To this purpose, we expand upon the stochastic multi-scale methodology developed in , where it was shown that such a system can be described by an age-structured birth-and-death process, instead of a branching process (Danesh et al., 2012;Durrett, 2013). The coupling between intracellular and the birth-and-death dynamics is carried out through a novel method to obtain the birth rate from the stochastic cell-cycle model, based on a mean-first passage time approach. Cell proliferation is assumed to be activated when one or more of the proteins involved in the cell-cycle regulatory pathway hit a threshold. This view allows us to calculate the birth rate as a function of the age of the cell and the extracellular oxygen in terms of the associated mean-first passage time (MFPT) problem (Redner, 2001). The present approach differs from that in  in that our treatment of the intracellular MFPT is done in terms of a large deviations approach, the so-called optimal path theory (Freidlin and Wentzell, 1998;Bressloff and Newby, 2014).
This methodology allows us to explore the effects of intrinsic fluctuations within the intracellular dynamics, in particular a model of the oxygen-regulated G1/S which dictates when cells are prepared to divide, as a source of heterogeneous behaviour: fluctuations induce variability in the birth rate within the population (even to the point of rendering some cells quiescent, i.e. stuck in G0) upon which a cell-cycle dependent therapy acts as a selective pressure.
This paper is organised as follows. Section 2 provides a summary of the structure of the multi-scale. In Section 3, we give a detailed discussion of the intracellular dynamics, i.e. the stochastic model of the oxygen-regulated G1/S transition, and its analysis. In Section 4, we summarise the formulation of the age-structured birth-and-death process, the numerical simulation technique, and the mean-field analysis of a homogeneous population. In Section 5, we discuss how noise within the intracellular dynamics induces heterogeneity in the population and analyse the stochastic dynamics of competition for a limited resource within such heterogeneous populations. In Section 6 we further study the effects of noise-induced heterogeneity on the emergence of drug resistance upon administration of a cell cycle-specific therapy. Finally, in Section 7 we summarise our results and discuss our conclusions as well as avenues for future research.

Summary of the multi-scale model
Before proceeding with a detailed discussion of the different elements involved in the formulation of the stochastic multi-scale model, it is useful to provide a general overview of the overall structure of the model, which is closely related to that of the model proposed in Alarcón et al. (2005).
The model we present in this article integrates phenomena characterised by different time scales, as schematically shown in Fig. 1. This model intends to tackle the growth and competition of cellular populations under the restriction of finite amount of available resources (in this case, oxygen) supplied at a finite rate,S.
The general approach used in this model is a natural generalisation of the standard continuous-time birth-and-death Markov process and its description via a Master Equation (Gardiner, 2009). As we will see, the consideration of the multi-scale character of the system, i.e. the inclusion of the physiological structure associated with the cell-cycle variables, introduce an age-structure within the population: the birth rate depends on the age of cell (i.e. time elapsed since last division) which determines, through the corresponding cell-cycle model, the cell-cycle status of the corresponding cells.
The evolution of the concentration of oxygen, c(t), (resource scale, see Fig. 1) is modelled by: where N T is the number of cellular types consuming the resource c, and N i (t), = … i N 1, , T , is the number of cells of type i at time t. Note that, in general, N i (t) is a stochastic process, and, therefore, in principle Eq. (1) should be treated as a stochastic differential equation (Oksendal, 2003).
The second sub-model considered in our multi-scale model, associated with the intracellular scale (see Fig. 1), is a stochastic model of oxygen-regulated cell-cycle progression. This sub-model is formulated using the standard techniques of chemical kinetics modelling (Gillespie, 1976) so that the mean-field limit of the stochastic model corresponds to the deterministic cell-cycle model formulated in Bedessem and Stephanou (2014). This model provides the physiological state of the cell in terms of the number of molecules of each protein involved in the cell-cycle of a cell of a given age, a. From such a physiological state, we derive whether the G1/S transition has occurred. The cell-cycle status of a cell of age a is determined in terms of whether the abundance of certain proteins which activate the cell-cycle (cyclins) have reached a certain threshold. In our particular case, if at age a, the cyclin levels are below the corresponding threshold, the cell is still in G1. If, on the contrary, the prescribed threshold level has been reached, the cell has passed onto S, and therefore is ready to divide. This implies that the probability of a cell having crossed the threshold of cyclin levels at age a can be formulated in terms of a mean firstpassage time problem (MFTP) in which one analyses the probability of a Markov process to hit a certain boundary (Redner, 2001;Gardiner, 2009). Unlike our approach in , based on approximating the full probability distribution of the stochastic cell cycle model, in the present approach, passage time is (approximately) solved in terms of an optimal trajectory path approach (Freidlin and Wentzell, 1998;Bressloff and Newby, 2014). At the interface between the intracellular and cellular scales sits our model of the (age-dependent) birth rate, which defines the probability of birth per unit time (cellular scale) in terms of the cell cycle variables (intracellular scale). The rate at which our cell-cycle model hits the cyclin activation threshold, i.e. the rate at which cells undergo G1/S transition, is taken as proportional to the birth rate. In particular, the birth rate is taken to be a function of the age of the cell as well as the concentration of oxygen, as the oxygen abundance regulates the G1/S transition age, ( ) a c G S 1/ , i.e. the time (age) elapsed between the birth of a cell and its G1/S transition: i.e. cell division occurs at a constant rate, τ − p 1 , provided cells undergo the G1/S transition and H is the Heaviside function. In other words, we consider that the duration of the G1 phase is regulated by the cell cycle model, whereas the duration of the S-G2-M is a random variable, exponentially distributed with average duration equal to τ p (see Fig. 1).
The third and last sub-model is that associated with the cellular scale. It corresponds to the dynamics of the cell population and is governed by the Master Equation for the probability density function of the number of cells (Gardiner, 2009). The stochastic process that describes the dynamics of the population of cells is an age-dependent birth-and-death process where the birth rate is given by Eq. (2) where ( ) a c G S 1/ is provided by the intracellular model. The death rate is, for simplicity, considered constant. As a consequence of the fact that the birth rate is age-dependent, our Multi-Scale Master Equation (MSME) does not present the standard form for unstructured populations, rather, it is an age-dependent Master Equation.
A detailed description of each sub-model and its analysis is given in Sections 3 and 4.
3. Intracellular scale: stochastic model of the oxygen-regulated G 1 /S transition

Biological background
Cell proliferation is orchestrated by a complex network of protein and gene expression regulation, the so-called cell cycle, which accounts for the timely coordination of proliferation with growth and, by means of signalling cues such as growth factors, tissue function (Yao, 2014). Dysregulation of such an orderly organisation of cell proliferation is one of the main contributors to the aberrant behaviour observed in tumours (Weinberg, 2007).
The cell cycle has the purpose of regulating the successive Schematic representation of the different elements that compose our multi-scale model. We show the different levels of biological organisation as well as associated characteristic time scales  associated to each of these layers: resource scale, i.e. oxygen which is supplied at a constant rate and consumed by the cell population, cellular scale, i.e. oxygen-regulated cell cycle progression which determines the age-dependent birth rate into the cellular layer, and, finally, the cellular scale, which is associated to the stochastic population dynamics.
activation of the so-called cyclin-dependent kinases (CDKs) which control the progression along the four phases of cycle: G1 (first gap phase), S (DNA replication), G2 (second gap phase), and M (mitosis) Goldbeter, 2009, 2011;Gerard et al., 2015). These four phases must be supplemented with a fifth, G0, which accounts for cells that are quiescent due to lack of stimulation (i.e. absence of growth factors, lack of basic nutrients, etc.) to proliferate. Recent models of the cell cycle organise the complex regulatory network into CDK modules, each centred around a cyclin-CDK complex which is key for the transition between the cell cycle phases (see, for example, Gerard and Goldbeter, 2009): cyclin D/CDK4-6 and cyclin E/CDK2 regulate progression during the G1 phase and elicit the G1/S transition, cyclin A/CDK2 promote progression during S phase and orchestrates the S/G2 transition, and, finally, cyclin B/CDK2 brings about the G2/M transition. The activity of each of these cyclin-CDK complexes is regulated in a timely manner, so that each phase of the cell cycle ensues at the proper time, by means of transcriptional regulation, post-transcriptional modifications (e.g. phosphorylation), and degradation (via ubiquitination) in which a large number of other components participate, including transcription factors, enzymes, ubiquitins, etc.
In the present paper, we propose a coarse-grained description of the cell cycle phases by lumping S, G2, and M into one phase, so that we consider a two-phase model G1 and S-G2-M, as shown in Fig. 1, Alarcón et al. (2004). In particular, we consider that cells can only divide once they have entered the S-G2-M phase at a constant rate. Entry in S-G2-M is regulated by a (stochastic) model of the G1/S transition which takes into account the regulation of the duration of the G1 phase by hypoxia (lack of oxygen).
The abundance of oxygen is known to be one of the factors that regulates the duration of the G 1 phase of the cell cycle. The issue of the regulation of the G 1 /S transition by the oxygen concentration has been the subject of several modelling studies (Alarcón et al., 2004;Bedessem and Stephanou, 2014). These models focus on the hypoxia-induced delay of the activation of the cyclins either through activation of cyclin inhibitors (Alarcón et al., 2004) or via up-regulation of the HIF-1α transcription factor (Bedessem and Stephanou, 2014). From the modelling point of view, both of them are mean-field models, thus neglecting fluctuations. In this section, we formulate a stochastic version of the model of Bedessem and Stephanou (2014), of which a schematic representation is shown in Fig. 2.
HIF-1 mediates adaptive responses to lack of oxygen (Semenza, 2013). HIF-1 is a heterodimer consisting of two sub-units: HIF-1α and HIF-1β. Whilst the latter is constitutively expressed, HIF-1α is O 2 -regulated. In the presence of adequate oxygen availability is negatively regulated by the von Hippel-Lindau (VHL) tumour suppressor protein, which allows HIF-1α for degradation. VHL loss-of-function mutations are common in many types of tumours, which allows for de-regulated HIF-1α degradation (Semenza, 2013). HIF-1 is involved in a number of cellular responses including switch from oxidative phosphorylation to glycolysis, activation of angiogenic pathways, and inhibition of cell cycle progression (Bristow and Hill, 2008;Semenza, 2013).

Mean-field model of the regulation of the G1/S transition by hypoxia
Bedessem and Stephanou formulate a model of the effect of hypoxia (mediated by HIF-1) on the timing of the G1/S transition (Bedessem and Stephanou, 2014). The involvement of HIF-1 in cell cycle regulation is complex and not completely understood. There is evidence that HIF-1 delays entry into S(-G2-M) phase by activating p21, a CDK inhibitor (Koshiji et al., 2004;Ortmann et al., 2014). HIF-1 up-regulation of p21 mediates indirect down-regulation of cyclin E (Goda, 2003). Further to HIF-1 regulation of the cell cycle, there exists a feedback regulation of cell cycle proteins of HIF-1. Hubbia et al. (2014) report that CDK1 and CDK2 physically and functionally interact with HIF-1α: CDK1 down-regulates lysosomal degradation of HIF-1α, thus consolidating its stability and promoting its transcriptional activity. On the contrary, CDK2 activates lysosomal degradation of HIF-1α and promotes G1/S transition. Bedessem and Stephanou (2014) do not take into account all these issues and, for simplicity, chose to focus on the well-documented effect of HIF-1 on cyclin D (Wen et al., 2010;Ortmann et al., 2014). Bedessem and Stephanou (2014) model formulation is based on the following assumptions: 1. The G1/S is modelled by a biological switch which emerges from the mutual inhibition between a cyclin (in this case, cyclin E) and a ubiquitin complex (SCF complex): The latter marks the former for degradation whereas cyclin E mediates inactivation of the SCF complex. This mutual inhibition gives rise to a bistable situation in which two stable fixed points coexist, each associated with the G1 phase (high SCF activity, low cyclin E concentration) and the S-G2-M phase (low SCF activity, high cyclin E concentration). Activation and inactivation of the SCF complex are assumed to follow Michaelis-Menten kinetics. 2. As in previous models (Tyson and Novak, 2001;Novak and Tyson, 2004;Alarcón et al., 2004), the G1/S transition is brought about by triggering a saddle-node bifurcation, whereby the G1 phase fixed point collides with the unstable saddle, driving the system into S-G2-M fixed point. Two factors drive the system through this bifurcation: cell growth (logistic increase of the cell mass, Tyson and Novak, 2001) and activation of the E2F transcription factor. In Bedessem and Stephanou (2014), both factors are taken to up-regulate the transcription of cyclin E. 3. Activation of E2F is modelled in terms of the E2F-Retinoblastoma protein (RBP) switch (Lee et al., 2010). Briefly, E2F is captured (bound) by unphosphorylated RBP. Upon phosphorylation, RBP releases E2F which activates transcription of G1/Stransition promoting cyclins, such as cyclin E (Alberts et al., 2002). Further, E2F can be in unphosphorylated (active) form and phosphorylated (inactive) form. Following Novak and Tyson (2004), Bedessem and Stephanou assume that fraction of active E2F and RBP-bound E2F are in adiabatic equilibrium with unphosphorylated RBP and total free E2F (Bedessem and

HIF−1a CycD
RbNP E2FA CycE SCF Growth Hypoxia Fig. 2. Schematic representation of the elements involved in the model of hypoxiaregulated G 1 /S transition proposed by Bedessem and Stephanou (2014). Within the framework of this model, the negative-feedback between CycE and SCF is the key modelling ingredient for the system to emulate the behaviour of a cell during the G1/S transition. The relative balance between CycE (which promotes the G1/S transition) and SCF (G1/S transition inhibitor) is regulated by growth and hypoxia, so that the timing of the transition depends upon these two factors.
These basic hypotheses are primarily based on previous models (Tyson and Novak, 2001;Novak and Tyson, 2004;Alarcón et al., 2004). The resulting pathway is schematically represented in Fig. 2.

Stochastic G1/S transition model
Based on the hypotheses summarised in Section 3.2, Bedessem and Stephanou (2014) formulated a mean-field model which is able to reproduce such behaviours as delay of the G1/S due lo lack of oxygen as well as hypoxia-induced quiescence. Here, we present a stochastic model (see Fig. 3 and Table 1), whose mean field limit is the model formulated in Bedessem and Stephanou (2014). We analyse this model using the stochastic quasi-steady state approximation we have developed in Alarcón (2014) and de la Cruz et al. (2015). The deterministic model formulated in Bedessem and Stephanou (2014), which we have briefly described in Section 3.2, is based on a series of reactions shown in Fig. 3, which include Michaelis-Menten kinetics for activation and inactivation of SCF complexes. Our stochastic model of the G1/S transition builds upon the stochastic (Markovian) description of the same set of reactions.
Our model is predicated on the stochastic dynamics of the state vector being described by a Markov jump process (Grimmett and Stirzaker, 2001;Kampen, 2007), whereby the state of the system, X, changes by an amount r i when the elementary reaction i occurs. The waiting time between Markovian events is exponentially distributed, the process is characterised by the associated transition rates, i.e.
( ( + Δ ) = + | ( ) = ) = ( )Δ + (Δ ) P X a a X r X a X W X a O a i i 2 . Using law of mass action (Gillespie, 1976) as our basic modelling framework, the transition rates of each elementary process are given in Table 1. Once we have determined the transition rates associated with each elementary reaction (or channel), the dynamics of the system is given by the Chemical Master Equation of a (non-structured) Markov Process, ( ) X a : is the probability of the state vector of the system to be X at age, i.e. the time reckoned from the last division, a. The transition rates, ( ) W X i , the vectors r i (whose components are the variation of the number of each chemical species upon occurrence of reaction i) are given in Table 1 and determine the dynamics of the system.
Even for moderately complex models, Eq.
(3) has no solution in closed form. Therefore, in order to study the properties of the system one must resort to numerical simulation (Monte Carlo) or asymptotic approximations. In the next section, we present an asymptotic analysis based on a recently developed form of stochastic quasi-steady state approximation.   Bedessem and Stephanou (2014). The reactions correspond (from top to bottom) to: hypoxia-inhibited synthesis and degradation of CycD, enzyme-catalysed, CycE-mediated inactivation of SCF, enzyme-catalysed activation of SCF, synthesis (regulated by growth and active E2F) and degradation (up-regulated by active SCF) of CycE, synthesis and degradation of Rb, and, last, synthesis and degradation of E2F. The negative feedback (mutual inhibition) between SCF and CycE mediates bistable behaviour in this model of the G1/S transition. Some of the transition rates associated to the reactions shown in are not constant but depend on the number of molecules of chemical species i present at a particular time, X i . For the definition of the quantities X i , we refer the reader to Table 1. Table 1 Reaction probability per unit time, The mass is assumed to grow following a logistic law: . Furthermore, Following Novak and Tyson (2004), we assume that at each time, the active E2F, The equilibrium between E2F-Rb complexes; free E2F and free Rb is given by Novak and Tyson (2004): Variable Description X 1 , X 8 Number of Cyclin D and Cyclin E molecules, respectively X 2 , X 5 Number of inactive and active SCF molecules, respectively X 3 , X 6 Number of SCF-activating and SCF-inactivating enzyme molecules, respectively X 4 , X 7 Number of enzyme-inactive SCF and enzyme-active SCF complexes, respectively X 9 , X 10 Number of free RbP and E2F molecules, respectively we have developed a stochastic version of the classical QSS approximation, the so-called semi-classical QSS approximation (SCQSSA) which, within the framework of the optimal path theory, allows us to tackle systems which exhibit separation of time scales, such as enzyme-catalysed reactions. Since these types of reaction feature prominently in our stochastic model of the hypoxia-regulated G1/S transition (see Fig. 3), we will use the SCQSSA to analyse the effects of intrinsic noise on the stochastic model of the hypoxiaregulated G1/S transition (as determined by the transition rates shown in Table 1). This approximation allows us to study noiseinduced phenomena which are relevant for the timing of the G1/S transition and, therefore, bear upon the population dynamics.
Following the SCQSSA methodology summarised in Appendix A, we derive the following set of equations which describe the optimal path associated with the stochastic G1/S transition model, Table 1: where a is the rescaled age variable = a kESt 7 (see Appendix B, Table B1). The oxygen dependencies enter the model though the [ ] H dependent parameter κ 2 (see Tables B.3-B.5). The reader is referred to Appendix A for a summary of the method and Appendix B for a detailed derivation of Eqs. (4)-(12) which hereafter is referred to as the semi-classical quasi-steady state approximation (SCQSSA) of the stochastic cell-cycle model and to Alarcón (2014) and de la Cruz et al. (2015) for a full account of the SCQSSA methodology. Note that there are other approximations which do not require the QSSA assumption, such as the one described in Kurtz (1978) and Ball et al. (2006). Furthermore, since we will be interested in the random effects associated to the enzyme-regulated dynamics of SCF, encapsulated in the parameters p 3 and p 6 , we have taken ( ) = p a 1 i for all ≠ i 3, 6. This corresponds to analysing the marginal distribution integrating out all the stochastic effects associated to all X i with ≠ i 3, 6.

Stochastic behaviour of the G1/S model
We now proceed to study the behaviour of the stochastic model of the oxygen-regulated G1/S transition. We pay special attention to those aspects in which we observe a departure of the stochastic system from the mean-field behaviour. In particular, we highlight the effects of modifying the relative abundance of SCFactivating and inactivating enzymes, including the ability of inducing oxygen-independent quiescence.
3.5.1. The relative abundance of the SCF activating and inactivating enzymes controls the timing of the G1/S transition In references Alarcón (2014) and de la Cruz et al. (2015), we have shown that, under SCQSSA conditions, the momenta p 3 and p 6 , i.e. the momenta coordinates associated with the SCF-activating and inactivating enzymes, respectively, are determined by the probability distribution of their initial (conserved) number. In particular, if we assume that the initial number of SCF-activating and inactivating enzyme molecules, E 1 and E 2 , is distributed over a population of cells following a Poisson distribution with parameter E, we have shown that (de la Cruz et al., 2015): . With this in mind, we can analyse the effect of changing the relative concentration of SCF-activating and inactivating enzymes on the timing of the G1/S transition. Our results are shown in Figs. 4 and 5. Fig. 4 illustrates that, for a fixed oxygen concentration, the G1/S transition is delayed by depriving the system of SCFactivating enzyme: as the ratio = p p e e / / 3 6 1 2 of SCF activating and deactivating enzyme increases, the G1/S transition takes longer to occur. Then, Fig. 5 shows that the G1/S transition age ( ) a c p p , , G S 1/ 3 6 decreases when p p / 3 6 increases. Furthermore, increasing the oxygen concentration c from c¼ 0.1 to c¼ 1 shifts the curve towards lower transition ages ( ) a c p p , , G S 1/ 3 6 . Note that this prediction is beyond the reach of the mean-field limit Bedessem and Stephanou (2014).

Induction of quiescence
In view of the results of Section 3.5.1, we have proceeded to a more thorough analysis of the effect of varying the ratio p p / 3 6 , which we recall that, within the SCQSSA, is equal to the ratio between the abundance of SCF-activating and inactivating enzymes, on the behaviour of the SCQSSA system Eqs. (4)-(12). In particular, we have investigated the bifurcation diagram of Eqs.
(4)-(12) with p p / 3 6 as the control parameter. Our results are shown in Fig. 6. We observe that, regardless of the value of m, there exists a range of values of the control parameter for which the saddle-node bifurcation, which gives rise to the G1/S transition, does not occur (i.e. only the G1-fixed point is stable). This result implies that depletion of SCF-activating enzyme, or, equivalently, over-expression of SCF-inactivating enzyme can stop cell-cycle progression by locking cells into the so-called G0 state, i.e. quiescence.
These results are confirmed by direct simulation of the stochastic cell-cycle model (Table 1) using Gillespie's stochastic simulation algorithm (Gillespie, 1976), see Fig. 7.
3.6. Scaling theory of the G1/S transition age We finish our analysis of the intracellular dynamics by formulating a scaling theory of one of the fundamental quantities in our multi-scale model, namely, the G1/S transition age, ( ) a c p p , , G S 1/ 6 3 , which determines the age-dependent birth rate (see Eq. (2)). In this section, we will show that, in spite of the complexity of the SCQSSA formulation of the oxygen-dependent cell-cycle progression model (see Eqs. (4)-(12)), a G S 1/ exhibits remarkable regularities with respect to its dependence on the oxygen concentration and the cell-cycle parameters p 6 and p 3 . Such regularities hugely simplify our multi-scale methodology.
We can see this regularity in Fig. 8. Fig. 8(a)  for small c, with the disagreement increasing with increasing c. Thus, a good scaling approximation for a G S 1/ is given by  Fig. 8(d). The critical oxygen concentration for quiescence c cr can be estimated analytically (with parameter values taken from Table B2):  Table B2. The parameter a 0 can be estimated as follows. Let A(c) be defined as: goes from 3 to 1. The parameters b i , e i , and J i are defined in Appendix B, Table B2 Fig. 4. Series of plots illustrating how the ratio p p / 6 3 , which is associated with the ratio of the number of SCF-inactivating and SCF-activating enzymes, modulates the timing of the G1/S transition. Parameter values as given in Table B2. Initial conditions are provided in Table B4.

Cellular scale: multi-scale Master Equation and path integral formulation
We start by summarising the formulation of the multi-scale Master Equation (MSME) for the population dynamics model in terms of an age-structured stochastic process . First, we consider a simple age-dependent birthand-death process where ( ) n a t , stands for the number of cells of age a at time t. Both a and t are dimensionless according to the scaling prescribed in Table B1, Appendix B. The time variable is → t kESt 7 and a is defined after Eq. (12). The offspring of such cells with age a¼ 0.
where ν ( ( ) ) = ( + ( )) ( ) W n a a t b a n a , , (see Table 2), μ is the (age-independent) death rate, and the age-dependent birth rate, b(a), is given by is the generalised coordinate associated with the concentration of CycE which must exceed a threshold value, CycE T for the cell-cycle to progress beyond the G 1 /S transition. Before this transition occurs, cells are not allowed to divide.
By re-arranging Eq. (17) and taking the limit δ δ = → t a 0, we obtain: where (see Table 2). Eqs. (19) needs to be supplemented with the appropriate boundary condition at a ¼0, ( = ) P n a t , 0, 0 . We proceed by first considering the number of births that occur within the age group = a a j during a time interval of length δt. Since we are assuming that our stochastic model is a Markov process where, within each age group, birth and death occur independently and with exponentially distributed waiting times, the number of births, ( ) B a j , is distributed according to a Poisson distribution: Its probability density is therefore given by: Since Eq. (20) is a convolution, using the well-known property of the probability generating function, the generating function associated with P(B), is given by Grimmett and Stirzaker (2001): , , is the generating function associated with ( ( )) P B a j j : j b a n a t p 1 j j Therefore, by taking δ δ . 21

Numerical method
Before proceeding further, we briefly describe the numerical methodology that we use to simulate the stochastic multi-scale model. In essence, the numerical method is an extension of the hybrid stochastic simulation algorithm used in  to accommodate the age structure of the cell populations we deal with here . For simplicity, we restrict our description of the algorithm to a homogeneous cell population. Its generalisation to heterogeneous populations composed of a variety of cellular types is straightforward.
Similarly to the procedure described in , we start by defining the population vector t is the set of indexes which label all those age groups which, at time, t, are represented within the population, i.e. all those age groups such that ( = ) > n a a t , 0 j . Having defined ( ) t , we summarise the numerical algorithm: 1. Set initial conditions: . We also set the value of the ratio of the SCF-regulating enzymes p p / 6 3 . 2. At some later time, t, the system is characterised by the quantities c(t), ( ) t , ( ) t , and N(t) 3. Generate two random numbers, z 1 and z 2 , uniformly distributed in the interval ( ) 0, 1 4. Use z 1 to calculate the exponentially distributed waiting time to the next event: 1 2 j j . The rates W 1 and W 2 are given in Table 2 5. Update the oxygen concentration at τ ( + ) c t by solving the associated ODE (1) between τ ( + ) t t , taking c(t) as initial condition. We use a 4 stage Runge-Kutta solver 6. Update the age-dependent birth rate for each ∈ j , i.e. Steps 3 to 10 are repeated until some stopping condition (e.g. ≥ t T) is satisfied Fig. 8. Series of plots showing the scaling analysis of the G1/S transition age. Plot (a) shows that below the critical value r cr , which corresponds to the value of the ratio p p / 6 3 , above which there is no transition to quiescence (i.e. if > p p r / cr 6 3 the G1/S transition age is finite when c¼ 0), ( ) a c p p , , G S 1/ 6 3 follows an algebraic decay with a universal p p / 6 3 -independent exponent provided that the oxygen, c, is rescaled by the critical oxygen concentration, 6 3 decays exponentially with the oxygen concentration with a characteristic concentration c 0 which, provided p p / 6 3 is larger enough than r , cr is p p / 6 3 -independent. Plot (b) shows how c cr varies as p p / 6 3 is changed. Similarly, plot (d) shows how + a varies p p / 6 3 changes. Parameter values as given in Table B2.

Table 2
This table summarises the elementary events involved in the age-dependent birthand-death process. Cells of age a produce offspring at a rate proportional to the age dependent birth rate, ( ) b a , Eq. (18). We are assuming that upon cell division both cells are reset to = a 0. Therefore, upon proliferation, one cell is removed from the population of age a and two cells are added to the population of age = a 0. For simplicity, death is assumed to be age-independent and to occur at a rate proportional to the death rate, μ.

Event
Reaction Transition rate,

Note that
Step 5 does not involve the use a stochastic method of integration of the ODE which rules the time evolution of the oxygen concentration, Eq. (1). This is due to the fact that, within the time interval τ ( + ] t t , , the population stays constant, so that Eq. (1) can be solved by means of a non-stochastic solver.

Steady-state of a homogeneous population: mean-field analysis
Before proceeding further, in order to check the numerical algorithm proposed in Section 4.1, we analyse how it compares with results regarding the steady-state of the mean-field limit of a homogeneous (i.e. composed by one cellular type only) population Hoppensteadt (1975). The mean-field equations associated with the stochastic multi-scale model are: with boundary condition: t baQ a td a 0, 2 , . and birth rate given by: 6 3 is given by Eq. (13) According to Hoppensteadt (1975), in order to ascertain whether a steady-state solution, i.e. whether the system settles onto an age distribution where the proportion of cells of each age does not change, we seek for a separable solution: . Defining μ ν ( ) = + ( ) a ba, we obtain: , is therefore given by: For this quantity to be positive τ ν < 1 p must hold. This condition states that for a steady state to be reached, the average waiting time to division after the G1/S transition, τ p , must be smaller than the average life span of the cell, ν À 1 . Eq. (29) determines the stationary value of the oxygen concentration, ∞ c . The long-time dynamics of Eq. (23) can therefore be summarised as follows: given the values of τ p , ν, and the cell-cycle parameters (see Table B2 in Appendix B), the population evolves and consumes oxygen until the oxygen concentration reaches a steady value ∞ c . At this point, the population of resident cells has settled onto a steady state where its age structure does not change. The total number of cells is also constant and given by (see Eq. (23)): In order to verify the age-structured SSA proposed in Section 4.1, we compare its results with the mean-field predictions, which should be in agreement with the stochastic behaviour of the system for large values of the carrying capacity, ∞ N . Results are shown in Fig. 9. We observe that, as predicted by our steady-state analysis, the stochastic simulations show how the resident population goes through an initial (oxygen-rich) phase of exponential growth. As the population grows, oxygen is depleted and the resident population eventually saturates onto a number of cells which fluctuates around the mean-field prediction of the carrying capacity (Eq. (30)).
Further verification of the validity of our numerical methodology is provided in Fig. 10. Our mean-field theory predicts that, everything else remaining unchanged, the average steady-state population should increase linearly with the rate of oxygen delivery, S (see Eq. 30). By contrast, the average equilibrium oxygen concentration does not depend on S, as shown in Eq. (29), and, consequently, it should stay constant upon increasing the rate of oxygen delivery. Fig. 10 shows that our simulations agree with these mean-field predictions.

Quasi-neutral competition within heterogeneous populations
A problem of fundamental importance in several biological and biomedical contexts is that of a population composed by a heterogeneous mixture of coexisting cellular types. A particularly relevant example of such a situation is that of cancer, where heterogeneity within the cancer cell population is assumed to be a major factor in the evolutionary dynamics of cancer as well as the emergence of drug resistance (Gatenby and Vincent, 2003;Merlo et al., 2006;Gillies et al., 2012;Greaves and Maley, 2012). Within this context, we are interested in (i) exploring the stability of a coexisting heterogeneous population and (ii) the effects such heterogeneity has on the long-term effects a cell-cycle dependent therapy.
The origins of heterogeneity of cancer cell populations is normally attributed to genetic variability arising from chromosomal instability and increased mutation rate (Merlo et al., 2006). Our discussion of the model of the G1/S transition (see Section 3) suggests a different source of heterogeneity associated with stochastic effects due to intrinsic fluctuations. We have shown, by means of both the SCQSSA analysis and stochastic simulations, that the relative abundance of SCF-regulating enzymes, which is quantified in our analysis by the ratio p p / 6 3 (see Sections 3.5.1, 3.6 and Appendix A). In Section 3.5.1 we have shown that the timing of the G1/S transition, and, consequently, the overall birth rate is strongly affected by changes in this quantity. In view, of this we associate heterogeneity to a distribution of the abundance of such enzymes, whereby we associate cell phenotypes with different values of the ratio p p / 6 3 . We further assume that this heterogeneity is hereditary, i.e. daughter cells inherit the value of the ratio p p / 6 3 from their mother.

Competition between two sub-populations
To proceed further, we consider the case of the birth-and-death dynamics of a heterogeneous population composed by two subpopulations, ( ) n a t ,  In the remaining part of this section we will consider two cellular populations, resident and invader. Each of these phenotypes are determined in terms of the values four parameters. The resident cells are characterised by two population-dynamics parameters, namely, the average time to division after the G1/S transition, τ p 1 , and the death rate, ν 1 . We consider two further parameters, p 3 1 and p 6 1 , associated with the cell-cycle progression dynamics of the resident cells (see Section 3.4). Similarly, the invader is characterised by the corresponding parameters: τ p 2 , ν 2 , p 3 2 and p 6 2 .

Mean-field coexistence versus quasi-neutral stochastic competition
We proceed to analyse the conditions under which two populations are capable of long-term coexistence. In particular, we analyse a scenario in which the mean-field description predicts long -term coexistence between within a heterogeneous population leads to mutual exclusion of all strands but one through socalled quasi-neutral stochastic competition (Lin et al., 2012;Kogan et al., 2014;. The starting point for our study is the mean-field analysis carried out in Section 4.2. According to these results, (mean-field) populations evolve until a concentration of oxygen, ∞ c , is reached so that the associated replication number The replication number of either population is given by:  This mean-field scenario is the basis for the study of the longterm stochastic dynamics of two populations which satisfy Eq. . The reproduction number is the average number of offspring per cell. We know from elementary considerations (Grimmett and Stirzaker, 2001) that the value of such quantity allows us to classify birth-death/branching processes. If > R 1 0 the population grows, on average, exponentially and has a finite probability of eventual survival. In this case, the process is referred to as super-critical. If < R 1 0 the population undergoes exponential decline (on average) and the extinction probability is equal to 1. The last case, in which = R 1 0 , the so-called critical case, the population, on average, stays constant. However, due to effects of noise, extinction occurs with probability 1, with the probability of survival up to time t asymptotically tends to ( ) ∼ − P t t S 1 (Kimmel and Axelrod, 2002). In our case, both the resident and the invader undergo a critical stochastic dynamics, where, once the steady state has established itself, the population evolves very close to the mean-field line of fixed points until fixation of one of the species (and, consequently, extinction of the other) occurs.
In order to numerically check this scenario, we could proceed to estimate the survival probability at time t, P S (t). However, since this quantity exhibits a fat-tail behaviour, this would be computationally costly. A more efficient method is to resort to the asymptotic of the extinction time with system size, which in this case can be identified with the carrying capacity, K (Hidalgo et al., 2015). Typically, a quasi-neutral competition is associated with an algebraic dependence of the average extinction time of either population, T E , on the system size, in this case determined by the carrying capacity, K (Doering et al., 2005;Lin et al., 2012;Kogan et al., 2014). In Fig. 11 we plot simulation results for the competition between two identical populations. In particular, we study how the average extinction time of either population, T E , varies as the carrying capacity is changed. We observe that this quantity exhibits a linear dependence on the carrying capacity: Our scenario produces the same qualitative results as Lin et al. (2012) and Kogan et al. (2014) who studied the average extinction time of birth-and-death processes engaged in quasi-neutral competition. In both papers, the average extinction time was reported to depend linearly on system size.

Study of the effects of cell-cycle-dependent therapy
In Section 5 we have analysed how the dependence of the cellcycle progression on the concentration of SCF-activating and SCFinactivating enzymes allows us to engineer heterogeneous populations where invasion and coexistence may occur. In this section, we further explore the ability of inducing quiescence by varying the ratio between SCF-activating and SCF-inactivating enzymes, this time in connection with the ability of such populations to withstand the effects of cell-cycle-dependent therapy (Powathil et al., 2012;Gabriel et al., 2012;Billy and Clairambault, 2013;Powathil et al., 2013).
In particular, we consider a scenario where two cellular populations coexist. Initially, one of these populations consists of a set of cells actively progressing through the cell-cycle which have reached a steady state characterised by the mean-field equilibrium Eqs. (28)-(30). The second population consists of quiescent cells whose cell-cycle is locked into the G0 phase and therefore do not proliferate. These cells are further assumed to undergo apoptosis at a very slow rate. More specifically, the ratio SCF-activating and SCF-inactivating enzymes in the quiescent cell population is such that, for the steady-state level of oxygen for the active cells (see Eq. (29)), cells are locked into the G1-fixed point (see Fig. 6).
In this section we show that the presence of a quiescent population within a heterogeneous population may lead to resistance to cell-cycle-dependent therapy. In particular, we show that, whereas such therapies effectively reduce or even eradicate the active cell population, the feedback between the therapy-induced decrease in cell numbers and the associated increase in oxygen availability can yield to the quiescent population to enter the active state and thus regrow the population. In this sense, we claim that the quiescent population has a stem-cell-like effect whereby, under the action of therapeutic agent, can repopulate the system Alarcón and Jensen, 2010.

Mean-field analysis
We start our mean field analysis by considering a heterogeneous population composed by cells of two types: type 1 and type 2 cells. Type 1 consists of cells with values of ≡ p p so that they are locked in G0 (i.e. not cycling).
The associated mean-field dynamics are given by: ( ) N t 2 are the total cell population of type 1 and type 2 cells, respectively: Fig. 11. Simulation results corresponding to the competition between the sub-populations of a heterogeneous populations. We show how the average extinction time of either population, T E , varies as the carrying capacity, K, changes. We observe that the dependence is linear. For simplicity, the two sub-populations are assumed to be identical (i.e. with the same characteristic parameter values for both the intracellular dynamics (cell-cycle) and the population-level dynamics (birth and death rates)). The carrying capacity = 4 , which correspond to = · · · · · K 0.9969 10 , 1.9301 10 , 2.6956 10 , 3.4367 10 , 4.2661 10 For simplicity, we assume that both populations consume oxygen at the same rate, k 1 . We further assume that ν ν ⪡ 2 1 , i.e. type 2 (quiescent) cells die at a much slower rate than type 1 (active) cells. In Section 4.2, we have already analysed under which conditions Eqs. which determines the steady-state value of the oxygen concentration ∞ c . At this point, the population of resident cells has settled onto a steady state where its age structure does not change. The total number of cells is also approximately constant and given by: where we have used the fact that ν ν ⪡ 2 1 . The cell-cycle parameters p 3 2 and p 6 2 have been chosen so that, for = ∞ c c , , i.e. type 2 cells are initially locked into G0 (see Fig. 6).
Once the population reaches this therapy-free quasi-equilibrium state, we assume that a therapy which only acts on proliferating cells is administered. Examples of such therapies abound in cancer treatment and can take the form of cell-cycle specific drugs or radiotherapy (Powathil et al., 2012;Gabriel et al., 2012;Billy and Clairambault, 2013;Powathil et al., 2013). We characterise the efficiency of the therapy by the so-called survival fraction, F S , i.e. the percentage of cells which survive the prescribed dose. For example, in radiotherapy F S is usually taken to be given by the linear quadratic model: D stands for the radiation dosage expressed in Grays and α and β are cell type-specific parameters. In the present context, we do not specify any particular form of therapy and we simply take ∈ [ ) F 0, 1 S . Initially, the therapy only affects type 1 cells (since type 2 are not proliferating). Therapy affects the birth and death rates of the type 1 population, which now read: The resulting mean-field equation is given by: The action of therapeutic agent on the active population initially induces a decline of the population, which, in turn, involves an increase in the available oxygen concentration. The latter has the effect of accelerating the rate of progression of type 1 cells through the cell cycle. In the absence of the quiescent population, eventually both effects would find a balance and the population of active cells would settle onto a new equilibrium characterised by:  since a G S 1/ is a decreasing function of the oxygen concentration. Reoxygenation during cell-cycle dependent therapy, in particular, radiotherapy, has been predicted by other models (Kempf et al., 2015).
Consider now the effect of this process on the type 2 cell population which is initially quiescent. We know that hypoxia-induced arrest of the cell cycle is reversible (Alarcón et al., 2004;Bedessem and Stephanou, 2014), i.e. upon increase of the concentration of oxygen quiescent cells may re-enter the cell cycle and become proliferating. Re-entry of quiescent cells into the cell cycle is predicated upon a sufficient increase in the oxygen contraction: > ( ) c c p p , , i.e the initial oxygen concentration is such that type 2 cells are quiescent, but, as the therapy is administered and proceeds to act upon the type 1 cells, the oxygen concentration increases until it reaches its critical re-entry concentration. At this point, type 2 cells abandon quiescence and become active and competition between type 1 and type 2 cells ensues. In order to assess the long-time behaviour of the system, we first study the equilibrium of the type 2 cell population upon re-entry into cell-cycle progression. Its mean-field dynamics is given by: and c is determined by Eq.
(38). Recall that, upon re-entering cell-cycle progression, type 2 cells are no longer immune to the therapy. Although in general the survival fraction is type-dependent, for simplicity we assume that F S has the same value for both cell types. The equilibrium condition is once again given in terms of the associated reproduction number, i.e. = R 1 0 T 2 which yields:

Critical dosage
In order to gain some degree of control over the behaviour described in the previous section, it would be useful to provide an estimate of the critical dosage above which therapy-induced reoxygenation is capable of activating quiescent cells. We characterise the therapy dose by means of the critical survival fraction, F S C . Recall that the characteristic equation for the oxygen-dependent growth rate, σ ( ) c 1 , of the population of active cells is given by: is satisfied, type 2 cells out-compete type 1 cells and the system, composed entirely of type 2 cells, settles onto the steady state prescribed by Eq. (52). In this sense, quiescent cells have stem-cell-like behaviour, in the sense that they can repopulate the system. In this second scenario, therapy is not completely without virtue since therapy drives the system to be taken over by a slower-cycling (less aggressive) phenotype.

Simulation results
In order to check the accuracy of the mean-field analysis carried out in Section 6.2 regarding the critical survival fraction for rescue from quiescence. We start by showing (Fig. 12) two typical realisations of the stochastic population dynamics which illustrate the rescue mechanism. In these simulations, we first let the active population settle on to its steady state. We then apply a sustained therapy with constant survival fraction. A more aggressive treatment (F S ¼0.6 in Fig. 12) greatly affects the active population: the amount of active cells killed by the therapy induces re-oxygenation of the population above the critical oxygen level for activation of the quiescent population whereupon the quiescent cells become proliferating. In Figs. 12(a) and (c), we show that upon activation of the quiescent population, a competition between both populations ensues, which eventually leads to extinction of the active population. A less aggressive therapy (F S ¼ 0.7 in Fig. 12) also induces death of the active population and re-oxygenation. However, in this case, the latter is not intense enough to induce activation of the quiescent cells (see Fig. 12(d)) and therefore the active cells will repopulate the system as the quiescent population stays on its course to eventual extinction, as shown in Fig. 12(b). Fig. 13 shows simulation results for the variation of probability of fixation of the quiescent population as the survival fraction of the therapy, F S , changes. Our simulation results show qualitative agreement with our mean-field theory (see Sections 6.1 and 6.2): as the survival fraction increases (i.e. the therapy becomes less efficient), the probability of fixation abruptly decreases from almost certainty of fixation to almost certainty of extinction. We observe that our mean-field theoretical predicts a critical value for F S slightly smaller than the observed when fluctuations due to finite size effects are present. However, we observe that, as the carrying capacity of the system is increased, the critical value of F S converges to the mean-field value.

Discussion
In this paper, we have presented and studied a stochastic multiscale model of a heterogeneous, resource-limited cell population. This model accounts for a stochastic intracellular dynamics (in this particular case, a model of the oxygen-regulated G1/S transition) and an age-structured birth-and-death process for the cell population dynamics. Both compartments are coupled by (i) a model for the time variation of resource (oxygen) abundance which regulates the rate of cell-cycle progression, and (ii) a model of the age-dependent birth rate which carries out the coupling between the intracellular and the cellular compartments (see Fig. 1 for a schematic representation of the model and Section 2 for a summary of the model formulation).
Our analysis of the stochastic dynamics of the oxygen-regulated G1/S transition, which is a generalisation of the mean-field model presented in Bedessem and Stephanou (2014), has revealed a number of previously unreported properties related to the presence of fluctuations. In particular, the optimal path theory and the quasi-steady state approximation allow us to explore the effect of the SCF-regulating enzymes on the timing of the G1/S transition. The relative abundance of SCF-activating and inhibiting enzymes regulates the rate at which cells reach the G1/S transition: excess of SCF-activating enzyme can delay the transition and even rendering the cell quiescent regardless of oxygen concentration beyond the predictions of the mean-field model (see Sections 3.5.1 and 3.5.2). Furthermore, we have shown that the effects on timing of the G1/S transition of the relative abundance of SCF-activating enzyme give rise to a scaling form dependence of the age to the transition, a G S 1/ , whereby this quantity is a function of ( ) a c p p , / G S 1/ 6 3 which takes the form of Eq. (13) (see Section 3.6). Taken together, these results imply that stochasticity in the intracellular dynamics naturally generates variability within the population: an otherwise homogeneous population presents a distribution of birth rates induced by variability in the relative abundance of SCF-activating enzyme within the population of cells.
This variation in the duration of the cell-cycle allows us to analyse the dynamics of a stochastic heterogeneous population under resource limitation conditions. To this end, we consider populations formed by sub-populations of cells characterised by differing relative abundance of SCF-activating and inhibiting enzymes. We further assume that this heterogeneity is heritable (i.e. daughter cells inherited the ratio of SCF-regulating enzymes from their mother). In this scenario, we have shown sub-populations within heterogeneous population engage in quasi-neutral competition: sub-populations of cells get extinct in an average time which is of the order of the carrying capacity of the system (see Section 5.2). In the context of modelling cancer cells populations, where heterogeneity is a main contributor to the complex dynamics of cancer (Anderson et al., 2006;Merlo et al., 2006;Anderson and Quaranta, 2008;Gillies et al., 2012;Greaves and Maley, 2012), this result is of relevance since it allows us to estimate the rate at which sub-populations or clones disappear from the tumour. The fact that this rate is proportional to the inverse of the carrying capacity reveals a highly dynamic scenario where clones are quickly decaying and being replaced within the tumour. This can have a profound impact on a variety of evolutionary phenomena such as emergence of drug resistant phenotypes. From the modelling perspective, we should note that quasi-neutral competition is a purely stochastic scenario: the mean-field limit predicts coexistence between the corresponding cell types.
We have further explored the issue of emergence of drug resistant cell types by analysing a case study in which a quiescent population can be rescued from latency by the application of a cell-cycle dependent therapy. Examples of such therapies are radiotherapy or cytotoxic drugs designed to target cells in specific stages of cell cycle progression (Powathil et al., 2012). In the particular example analysed in Section 6, we have shown that in mixed population composed by active (proliferating) and quiescent cells, if the drug is not efficient enough (characterised in terms of its associated survival fraction, F S ), quiescent cells are rescued from latency and eventually reach fixation within the population, i.e. the activated quiescent cells out-compete the original active cells until the latter population becomes extinct (Alarcón and Jensen, 2010). It is noteworthy the fact that in this process of quiescence rescue the original cancer population is replaced by a much more resistant population as the activated Our results show that the methods and models presented in this paper are of great potential importance for the analysis of the complex dynamics of heterogeneous populations under resource limitation, in particular for the study of emergence of drug resistance in heterogeneous cancer cell populations. Several important issues have been left out of the present work. A major contributor to heterogeneity within a cancer cell population is spatial heterogeneity (Anderson et al., 2006;Gillies et al., 2012) which is closely related to micro-environmental heterogeneity. A further issue which should be analysed in depth concerns the scaling properties of the age to the G1/S transition (see Section 3.6), in particular whether this is a general property of the cell-cycle dynamics or rather a specific attribute of the stochastic model presented here. A thorough analysis of these issues falls beyond the scope of the present paper and are left for future research.  cón, 2014;de la Cruz et al., 2015). In order to proceed further, we assume, as per the Briggs-Haldane treatment of the Michealis-Menten model for enzyme kinetics (Briggs and Haldane, 1925;Keener and Sneyd, 1998), that the species involved in the system under scrutiny are divided into two groups according to their characteristic scales. More specifically, we have a subset of chemical species whose numbers, X i , scale as: , whilst the remaining species are such that their numbers, X j , scale as: . Key to our approach is the fact that S and E must be such that: We further assume (de la Cruz et al., 2015) that the generalised coordinates, Q i , scale in the same fashion as the corresponding variable X i , i.e.  Dimensionless variables and parameters corresponding.

Table B4
Initial conditions used in simulation system (4)