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.

Introduction

Global cell traits and behaviour in response to stimuli, the so-called phenotype, results from a complex network of interactions between genes and gene products which ultimately regulates gene expression. These networks of gene regulation constitute non-linear, high-dimensional dynamical systems whose structure has been shaped up by evolution by natural selection, so that they exhibit properties such as robustness (i.e. resilience of the phenotype against genetic alterations) and canalisation (i.e. the ability 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.

To tackle such complexity, multi-scale models of biological systems as diverse as cardiac systems (Smith et al., 2004, McCulloch, 2009, Hand and Griffith, 2010, Land et al., 2013), systems of developmental biology (Schnell et al., 2008, Oates et al., 2009, Hester et al., 2011, Setty, 2012, Walpole et al., 2013), and tumour growth systems (Alarcón et al., 2005, Jiang et al., 2005, Ribba et al., 2006, Macklin et al., 2009, Owen et al., 2009, Preziosi and Tosin, 2009, Tracqui, 2009, Byrne, 2010, Lowengrub et al., 2010, Osborne et al., 2010, Rejniak and Anderson,, Deisboeck et al., 2011, Perfahl et al., 2011, Travasso et al., 2011, Durrett, 2013, Powathil et al., 2013, Szabo and Merks, 2013, Chisholm et al., 2015, Curtius et al., 2015, Scott et al., 2016, Jagiella et al., 2016) have been developed. 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, Spill et al., 2016).

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 and Zeisberg, 2006, Grivennikov et al., 2010), or the extracellular matrix and its physical properties (Spill et al., 2016), strongly influence cancer cell behaviour. Note that hypoxia is also known to change the tumour microenvironment (Spill et al., 2016). 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 Guerrero and Alarcón (2015), 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 Guerrero and Alarcón (2015) 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, .

Fig. 1

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 (Guerrero and Alarcón, 2015) 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.

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 is the number of cellular types consuming the resource c, and N(t), , is the number of cells of type i at time t. Note that, in general, N(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 first-passage 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 Guerrero and Alarcón (2015), 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, , 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, , 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 τ (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 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 3, 4.

Intracellular scale: stochastic model of the oxygen-regulated G1/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 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) (Gerard and Goldbeter, 2009, Gerard and Goldbeter, 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 G1 phase of the cell cycle. The issue of the regulation of the G1/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.

Fig. 2

Schematic representation of the elements involved in the model of hypoxia-regulated G1/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.

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 O2-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:

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.

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.

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/S-transition 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 Stephanou, 2014). Furthermore, (Bedessem and Stephanou (2014)) takes cyclin D to phosphorylate RBP.

Last, Bedessem and Stephanou (2014) assume that HIF-1α inhibits synthesis of cyclin D. Following experimental evidence reported in Wen et al. (2010), they assume that the level of HIF-1α depends exponentially of the oxygen concentration.

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, de la Cruz et al., 2015.

Fig. 3

Reactions for the stochastic version of the model proposed by 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. For the definition of the quantities X, we refer the reader to Table 1.

Table 1

Reaction probability per unit time, . The mass is assumed to grow following a logistic law: , where and a is the age of the cell (i.e. the time elapsed since birth). According to Bedessem and Stephanou (2014), the level of active HIF-1 , , depends exponentially on the extracellular oxygen concentration, : . Furthermore, Following Novak and Tyson (2004), we assume that at each time, the active E2F, , is the fraction of unphosphorylated free E2F factor, . The equilibrium between E2F–Rb complexes; free E2F and free Rb is given by Novak and Tyson (2004): .

Variable	Description
X1, X8	Number of Cyclin D and Cyclin E molecules, respectively
X2, X5	Number of inactive and active SCF molecules, respectively
X3, X6	Number of SCF-activating and SCF-inactivating enzyme molecules, respectively
X4, X7	Number of enzyme-inactive SCF and enzyme-active SCF complexes, respectively
X9, X10	Number of free RbP and E2F molecules, respectively

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 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. . 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, :where 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, , the vectors r (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.

Semi-classical quasi-steady state analysis of the stochastic G1/S transition model

In Alarcón, 2014, de la Cruz et al., 2015, Sanchez-Taltavull,, 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 hypoxia-regulated G1/S transition (as determined by the transition rates shown in Table 1). This approximation allows us to study noise-induced 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 (see Appendix B, Table B1). The oxygen dependencies enter the model though the dependent parameter κ2 (see Table B1, Table B2, Table B3). The reader is referred to Appendix A for a summary of the method and Appendix B for a detailed derivation of Eqs. (4), (B.36), (B.37), (7), (8), (9), (10), (11), (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, 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 Ball et al., 2006, Kurtz, 1978. Furthermore, since we will be interested in the random effects associated to the enzyme-regulated dynamics of SCF, encapsulated in the parameters p3 and p6, we have taken for all . This corresponds to analysing the marginal distribution integrating out all the stochastic effects associated to all X with .

Table B1

Dimensionless variables and parameters corresponding.

Dimensionless variables	Dimensionless parameters
a=k_7ESt	ϵ=E/S
q_1=Q_1/S	[e2f]^_tot=[e2f]_tot/S
q_2=Q_2/S	κ_1=k_1/(k_7ES^2)
q_3=Q_3/E	κ_2=k_2/(k_7ES^2)
q_4=Q_4/E	κ_3=k_3/(k_7ES)
q_5=Q_5/S	κ_4=k_4/(k_7S)
q_6=Q_6/E	κ_5=k_5/(k_7S^2)
q_7=Q_7/E	κ_6=k_6/(k_7S^2)
q_8=Q_8/S	κ_8=k_8/(k_7S)
q_9=Q_9/S	κ_9=k_9/(k_7S)
q_10=Q_10/S	κ_10=k_10/(k_7ES)
	κ_11=k_11/(k_7ES)
	κ_12=k_12/(k_7E)
	κ_13=k_13/(k_7ES^2)
	κ_14=k_14/(k_7ES)
	κ_15=k_15/(k_7E)
	κ_16=k_16/(k_7ES^2)
	κ_17=k_17/(k_7ES)

Table B2

Equivalence between the parameters of the stochastic model of the G1/S transition (see Table 1) and the parameters of the associated mean-field model as formulated by Bedessem and Stephanou (2014).

Parameters	Parameters
k_1=a_1*S	k_11=b_2
k_2=a_3*S*[H]	k_12=b_3/S
k_3=a_2	k_13=d_2*S
k_5=J_2*k_4S−k_6	k_14=d_2
k_6=e_1S/E	k_15=d_1/S
k_8=J_1*S−k_9	k_16=g_1[E2F]_totS
k_9=e_2/E	k_17=g_1
k_10=b_1

Table B3

Parameters values of the mean-field model of the hypoxia regulated G1/S transition proposed by Bedessem and Stephanou (2014).

Parameter	Value	Reference
a_1	0.51	Bedessem and Stephanou (2014)
a_2	1	Bedessem and Stephanou (2014)
a_3H_0	0.0085	Bedessem and Stephanou (2014)
b_1	0.018	Bedessem and Stephanou (2014)
b_2	0.5	Bedessem and Stephanou (2014)
b_3	1	Bedessem and Stephanou (2014)
d_1	0.2	Bedessem and Stephanou (2014)
d_2	0.1	Bedessem and Stephanou (2014)
e_1	1	Bedessem and Stephanou (2014)
e_2	14	Bedessem and Stephanou (2014)
m_⋆	10	Tyson and Novak (2001)
J_3,J_4	0.04	Tyson and Novak (2001)
g_1	0.016	Bedessem and Stephanou (2014)
[E2F]_tot	1	Bedessem and Stephanou (2014)
S	10
E	1
r_cr	0.04

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 SCF-activating and inactivating enzymes, including the ability of inducing oxygen-independent quiescence.

The relative abundance of the SCF activating and inactivating enzymes controls the timing of the G1/S transition

In references Alarcón, 2014, de la Cruz et al., 2015, we have shown that, under SCQSSA conditions, the momenta p3 and p6, 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, E1 and E2, 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 Fig. 4, Fig. 5. Fig. 4 illustrates that, for a fixed oxygen concentration, the G1/S transition is delayed by depriving the system of SCF-activating enzyme: as the ratio 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 decreases when increases. Furthermore, increasing the oxygen concentration c from c=0.1 to c=1 shifts the curve towards lower transition ages . Note that this prediction is beyond the reach of the mean-field limit Bedessem and Stephanou (2014).

Fig. 4

Series of plots illustrating how the ratio , 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.

Fig. 5

Plot showing how the G1/S transition age, , changes as the ratio , which is determined by the ratio of the (conserved) amounts of SCF activating and inactivating enzymes (Alarcón, 2014, de la Cruz et al., 2015), varies. We show for c=1 (blue circles) and c=0.1 (red squares). Parameter values as given in Table B2. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

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 , 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), (B.36), (B.37), (7), (8), (9), (10), (11), (12). In particular, we have investigated the bifurcation diagram of Eqs. (4), (B.36), (B.37), (7), (8), (9), (10), (11), (12) with 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.

Fig. 6

This figure shows the bifurcation diagram of the SCQSSA of the stochastic cell-cycle model for different values of the parameter m. The ratio is the control parameter. The order parameter is the steady state value of the generalised coordinate associated with active SCF, q5. Solid line corresponds to m=10, dash lines to m=8 and dotted lines to m=6. Parameter values as give in Table B3, , p=1 and c=1. Blue lines indicate stable steady state and red lines indicate unstable steady state. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

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.

Fig. 7

Simulation results for the stochastic model of the oxygen-regulated G1/S transition defined by the transition rates given in Table 1. We have plotted the probability , where T=100, with different values of . 1000 realisations and m=5.0.

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, , 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), (B.36), (B.37), (7), (8), (9), (10), (11), (12)), exhibits remarkable regularities with respect to its dependence on the oxygen concentration and the cell-cycle parameters p6 and p3. Such regularities hugely simplify our multi-scale methodology.

We can see this regularity in Fig. 8. Fig. 8(a) shows as a function of the ratio , for six different values of , where is the critical value. We see that all graphs fall together on a straight line in log–log space, indicating a power-law dependence of on . On the other hand, Fig. 8(c) shows the dependence of the normalized function , on c, for five values . Here, . We see that there is good agreement between the five values for small c, with the disagreement increasing with increasing c. Thus, a good scaling approximation for is given byHere, c0, , and β are constants. According to our analysis of the data presented in Fig. 8, these constants are given by , and . Likewise, is obtained by fitting to the data presented in Fig. 8(d). The critical oxygen concentration for quiescence c can be estimated analytically (with parameter values taken from Table B2):where a1, a3, d1, d2, , and H0 are parameters defined in Appendix B, Table B2. The parameter a0 can be estimated as follows. Let A(c) be defined as:where . , where c is the critical value of the oxygen concentration for which the saddle node bifurcation occurs, i.e. the critical value for which the number of real, positive solutions of the equation:goes from 3 to 1. The parameters b, e, and J are defined in Appendix B, Table B2

Fig. 8

Series of plots showing the scaling analysis of the G1/S transition age. Plot (a) shows that below the critical value r, which corresponds to the value of the ratio , above which there is no transition to quiescence (i.e. if the G1/S transition age is finite when c=0), follows an algebraic decay with a universal -independent exponent provided that the oxygen, c, is rescaled by the critical oxygen concentration, . Plot (c) shows that if, by contrast, decays exponentially with the oxygen concentration with a characteristic concentration c0 which, provided is larger enough than is -independent. Plot (b) shows how c varies as is changed. Similarly, plot (d) shows how varies changes. Parameter values as given in Table B2.

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 (Guerrero and Alarcón, 2015). First, we consider a simple age-dependent birth-and-death process where 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 and a is defined after Eq. (12). The offspring of such cells with age a=0.where (see Table 2), μ is the (age-independent) death rate, and the age-dependent birth rate, b(a), is given bywhere is the generalised coordinate associated with the concentration of CycE which must exceed a threshold value, CycE for the cell-cycle to progress beyond the G1/S transition. Before this transition occurs, cells are not allowed to divide.

Table 2

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

Event	Reaction	Transition rate, W_k(n,a,t)	r_k
Birth	n(a)→n(a)−1∅→n(a=0)=2	W_1(n,a,t)=b(a)n(a)	r_1=−1
Death	n(a)→n(a)−1	W_2(n,a,t)=μn(a)	r_2=−1

By re-arranging Eq. (17) and taking the limit , we obtain:where (see Table 2). Eqs. (19) needs to be supplemented with the appropriate boundary condition at a=0, . We proceed by first considering the number of births that occur within the age group during a time interval of length . 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, , is distributed according to a Poisson distribution:

The total number of births delivered by the whole population during , B, is . 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):where is the generating function associated with :Therefore, by taking Eq. (21) determines since and

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 Guerrero et al. (2015) to accommodate the age structure of the cell populations we deal with here (Guerrero and Alarcón, 2015). 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 Guerrero and Alarcón (2015), we start by defining the population vector , where 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 . Having defined , we summarise the numerical algorithm:

Set initial conditions: , , . We also set the value of the ratio of the SCF-regulating enzymes .

At some later time, t, the system is characterised by the quantities c(t), , , and N(t)

Generate two random numbers, z1 and z2, uniformly distributed in the interval

Use z1 to calculate the exponentially distributed waiting time to the next event: where . The rates W1 and W2 are given in Table 2

Update the oxygen concentration at by solving the associated ODE (1) between taking c(t) as initial condition. We use a 4 stage Runge–Kutta solver

Update the age-dependent birth rate for each , i.e. . The quantity is given by Eq. (13)

Use z2 to determine which event occurs (i.e. birth or death within the jth sub-population) at time : event l occurs with probabilitywhere and and

If the randomly chosen event is (i.e. cell proliferation within age group j), then , for all , and

If the randomly chosen event is (i.e. cell death within age group j), then , for all

Finally, we update the set , i.e. the set of age groups for which

Steps 3 to 10 are repeated until some stopping condition (e.g. ) is satisfied

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 , 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:and birth rate given by:with 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 , we obtain:with cnt. to be determined. is therefore given by:The value of the parameter σ is obtained by means of the characteristic equation obtained by introducing Eq. (25) with a=0 into the boundary condition:

After some algebra, the characteristic equation (26) reads:From Eq. (27) we obtain the condition for the system to be in equilibrium, i.e. . Substituting in Eq. (26):where R0 is the average number of offspring per cell at equilibrium: if the system grows exponentially, the system dwindles, and if the population remains constant. The equilibrium condition allows us to find the value of for which such equilibrium exists:For this quantity to be positive must hold. This condition states that for a steady state to be reached, the average waiting time to division after the G1/S transition, τ, must be smaller than the average life span of the cell, ν−1. Eq. (29) determines the stationary value of the oxygen concentration, .

The long-time dynamics of Eq. (23) can therefore be summarised as follows: given the values of τ, ν, 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 . 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, . 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)).

Fig. 9

Plots showing simulation results of the stochastic multi-scale dynamics of a cell population. These plots show how, in agreement with our steady-state analysis, the population evolves until it reaches a steady-state where the population of resident cells fluctuates around its associated carrying capacity Eq. (30). Colour code: blue lines show the total resident cell population at time t, N(t) (panels (a) and (b)). Green lines (panels (c) and (d)) show the associated oxygen concentration, c(t). The results shown in this figure correspond to a single realisation of the process. Parameter values: , , , in panels (a) and (c), and in panels (b) and (d). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

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.

Fig. 10

Plots showing simulation results of the stochastic multi-scale dynamics of a cell population. These plots show how, in agreement with our steady-state analysis, the average population at equilibrium increases linearly as the rate of oxygen delivery, S, is changed (plot (a)). Also, in accordance with our mean-field theory, the average, steady-state oxygen concentration remains unchanged as S increases. Parameter values: , , in panels (a) and (c), and in panels (b) and (d). Average has been performed over 100 realisations.

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 co-existing 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 (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 . We further assume that this heterogeneity is hereditary, i.e. daughter cells inherit the value of the ratio 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 sub-populations, and , competing by a common resource, c(t), which regulates the rate of progression of the cell-cycle of each cell type. The stochastic dynamics of the whole population is determined by the associated multi-scale master equation:where and , ν1 and are the (age-independent) death rate and the birth rate of the resident population, and ν2 and are the (age-independent) death rate and the birth rate of the invader. The quantities and are determined by the (oxygen-dependent) rate of cell-cycle progression of each sub-population: is given by Eq. (13). The concentration of resource (oxygen), c(t), is determined by the following ODE:where for i=1,2. The associated initial conditions are given by:

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, , and the death rate, ν1. We consider two further parameters, and , associated with the cell-cycle progression dynamics of the resident cells (see Section 3.4). Similarly, the invader is characterised by the corresponding parameters: , ν2, and .

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 so-called quasi-neutral stochastic competition (Lin et al., 2012, Kogan et al., 2014, Guerrero et al., 2015).

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, , is reached so that the associated replication number , i=1,2. The replication number of either population is given by:Therefore, our theory predicts that, provided that there exists such that:is satisfied, long-term coexistence ensues, since the whole system (oxygen, resident and invader) is able to evolve to a state where both populations are in equilibrium (in the sense that for both populations) with the same concentration of oxygen. Furthermore, the number of resident cells, N1, and the number of invaders, N2, satisfy:Thus, the mean-field theory predicts that there exist a continuous of fixed points. The eventual convergence on to a particular point along the line of fixed points Eq. (36) depends on the initial conditions. Eq. (36) can be further simplified by assuming , in which case , where is the carrying capacity.

This mean-field scenario is the basis for the study of the long-term stochastic dynamics of two populations which satisfy Eq. (35), which are equivalent to . 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 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 the population undergoes exponential decline (on average) and the extinction probability is equal to 1. The last case, in which , 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 (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(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, on the system size, in this case determined by the carrying capacity, K Doering et al., 2005, Kogan et al., 2014, Lin et al., 2012. 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, varies as the carrying capacity is changed. We observe that this quantity exhibits a linear dependence on the carrying capacity:

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, 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 by varying the death rates of both populations. The values of the death rates are , which correspond to , respectively. Averages are done over 500 realisations of the hybrid stochastic model.

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 cell-cycle progression on the concentration of SCF-activating and SCF-inactivating 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), (29), (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 and (see Section 3.4) such that cells are actively progressing through the cell-cycle. Type 2 cells are characterised by values of and so that they are locked in G0 (i.e. not cycling). The associated mean-field dynamics are given by:with boundary conditions: and are the total cell population of type 1 and type 2 cells, respectively:For simplicity, we assume that both populations consume oxygen at the same rate, k1. We further assume that , 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. (38), (39), (40) reach a steady state:where is the average number of offspring per cell of type 1 at equilibrium. The associated equilibrium value of is then given by:which determines the steady-state value of the oxygen concentration . 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 . The cell-cycle parameters and have been chosen so that, for , , 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, i.e. the percentage of cells which survive the prescribed dose. For example, in radiotherapy F is usually taken to be given by the linear quadratic model: where 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 . 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:or, equivalently:where is the equilibrium oxygen concentration for the type 2 population with therapy. Note that and therefore since is a decreasing function of the oxygen concentration. Re-oxygenation 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: , where the critical oxygen concentration for type 2 cells, , depends on the momenta and or, equivalently, on the concentration of SCF-activating and SCF-inactivating enzymes. Taking this property into account, one can devise a scenario in which , 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:where 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 has the same value for both cell types. The equilibrium condition is once again given in terms of the associated reproduction number, i.e. which yields:

Since , we have that (see Eqs. (49), (52)). The latter inequality implies that the equilibrium oxygen concentration for type 2 cell, , is such that . It is easy to argue that, in these conditions, the type 2 population out-competes the type 1 cells: for , the growth rate of type 1 cells is positive whereas the growth rate of the type 2 cell population is negative (see Eq. (27)). This implies that, upon application of therapy and provided that is satisfied, the type 1 population declines and it is replaced by the type 2 population.

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 re-oxygenation is capable of activating quiescent cells. We characterise the therapy dose by means of the critical survival fraction, . Recall that the characteristic equation for the oxygen-dependent growth rate, , of the population of active cells is given by:In order to activate the quiescent population the oxygen concentration must raise above the critical value . For the oxygen concentration to grow above this threshold so that the active cell population continues to decline thus allowing the oxygen concentration to keep on raising. Therefore the critical value is such that , i.e.If , the decrease in the active cell population is enough to provide enough oxygen for the quiescent population to become active.

This analysis implies that in heterogeneous populations which include quiescent sub-populations, the effect of cell-cycle-dependent therapy does not eradicate the population. Rather the following two scenarios are possible. If , type 2 cells do not become activated and become eventually extinct, and the type 1 cells settle onto the steady state prescribed by Eq. (49). If, by contrast, 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=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=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. 12

Stochastic simulation results showing typical realisations associated with rescue of quiescence cells (plots (a) and (c)) and recovery of the proliferating population (plots (b) and (d)) upon application of a cell-cycle dependent therapy. The efficiency of the therapy is characterised by the survival fraction, F. Quiescence rescue is achieved when the survival fraction is set to a value which falls below the critical threshold, Eq. (54). If (plots (a) & (c)), the cell killing triggered by the therapy is enough to re-oxygenate the population above the activation threshold of the quiescent cells. By contrast, if (plots (b) and (d)), re-oxygenation is not enough to rescue latent cells from quiescence. Parameter values: , , , , . The subindex “1” corresponds to the active population whilst the subindex “2” denotes quantities associated with the quiescent population. The critical oxygen (as defined in 3.5.2, 3.6) is for the active cells and for the quiescent cells. Colour code: I all of the panels in this figure, blue (red) lines correspond to the time evolution of the total number of proliferating (quiescent) cells and green lines, to the time evolution of the oxygen concentration. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Fig. 13 shows simulation results for the variation of probability of fixation of the quiescent population as the survival fraction of the therapy, F, changes. Our simulation results show qualitative agreement with our mean-field theory (see 6.1, 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 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 converges to the mean-field value.

Fig. 13

This figure shows stochastic simulation results regarding the variation of the probability of fixation of the quiescence population varies the efficiency of a cell-cycle dependent therapy changes. The efficiency of the therapy is measured in terms of the survival fraction F. The orange vertical line represents the mean-field, theoretical critical survival fraction (Eq. 54). We observe that as the carrying capacity of the system increases, which in these simulations is achieved by increasing the rate of oxygen supply, S, the results of the stochastic simulations tend towards the mean-field prediction. Parameter values: , , , , . The subindex “1” corresponds to the active population whilst the subindex “2” denotes quantities associated with the quiescent population. The critical oxygen (as defined in 3.5.2, 3.6) is for the active cells and for the quiescent cells. Each data point corresponds to the average over 500 realisations of the stochastic population dynamics. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)

Discussion

In this paper, we have presented and studied a stochastic multi-scale 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 3.5.1, 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, , whereby this quantity is a function of 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), 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 quiescent cells are less sensitive to the therapy than the original active cells.

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.

Table B4

Initial conditions used in simulation system (4), (B.36), (B.37), (7), (8), (9), (10), (11), (12) Fig. 4.

Parameter	Initial condition	Reference
CycD	0.1	Bedessem and Stephanou (2014)
SCF	0.9	Bedessem and Stephanou (2014)
Rb	1.0	Bedessem and Stephanou (2014)
E2F	0.1	Bedessem and Stephanou (2014)
m	5.0	Bedessem and Stephanou (2014)