Stochastic Multiscale Models of Cell Population Dynamics: Asymptotic and Numerical Methods

. In this paper we present a new methodology that allows us to formulate and analyse stochastic multiscale models of the dynamics of cell populations. In the spirit of existing hybrid multiscale models, we set up our model in a hierarchical way according to the characteristic time scales involved, where the stochastic population dynamics is governed by the birth and death rates as prescribed by the corresponding intracellular pathways (e.g. stochastic cell-cycle model). The feed-back loop is closed by the coupling between the dynamics of the population and the intracellular dynamics via the concentration of oxygen: Cells consume oxygen which, in turn, regulate the rate at which cells proceed through their cell-cycle. The coupling between intracellular and population dynamics is carried out through a novel method to obtain the birth rate from the stochastic cell-cycle model, based on a mean-ﬁrst 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 corresponding mean-ﬁrst passage time. We then proceed to formulate the stochastic dynamics of the population of cells in terms of an age-structured Master Equation. Further, we have developed generalisations of asymptotic (WKB) methods for our age-structured Master Equation as well as a τ − leap method to simulate the evolution of our age-structured population. Finally, we illustrate this general methodology with a particular example of a cell population where progression through the cell-cycle is regulated by the availability of oxygen.

The interest driving the increased effort devoted to the development of multiscale models and techniques is motivated by the realisation that the magic bullet approach [53] to addressing complex diseases, in particular, in the case of cancer, might not be as effective as it once was thought.Regarding cancer treatment, this concept was put forward by Paul Ehrlich [53] and consists of a targeted therapy which acts specifically on cancer cells, leaving normal cells unharmed.With the arrival of the genomic area, this approach was expected to be considerably boosted.In fact, for over thirty years, research in oncology has been dominated by a genocentric approach where targeted therapies, i.e. drugs developed to interfere with specific cancer gene products have been the focus and ultimate aim of cancer biology [39].Advances in genomics and other omics (proteomics, epigenomics, etc.) have further boosted this approach.However, the success of this approach in terms of the development of new, efficient cancer drugs has felt short of expectations [39].
There are a number of reasons why the magic bullet approach has had limited success.Global cell traits and behaviour in response to stimuli, i.e. the phenotype, emerge from a complex network of interactions between genes and gene products which ultimately regulates gene expression (see, for example, the recent work by Lignet et al. regarding the VEGF signalling network [32]).These networks of gene regulation constitute non-linear, high-dimensional dynamical 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 [29].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 [2,14,34,41,44,48].These and other factors contribute towards a highly dynamics in biological tissues.In particular, because of all the layers of complexity involved, it is very difficult to assess the main tenet of the magic bullet approach, i.e. whether a therapeutic agent is going to be effective against the tumour and harmless to the surrounding healthy tissue.
In order to try and address these issues, a great deal of research has been done in the development and analysis of multiscale models.Such models are capable of incorporating within a single model different sub-models corresponding to different levels of biological organisation (intracellular, cell-to-cell interaction, whole-tissue level, etc.), which are usually characterised by diverse time and length scales, and the coupling between them, so that the global tissue behaviour can be analysed as an emergent property of the different coupled elements [11,14,33,41,47,56].
Multi-scale models can be formulated in a number of different ways.One of these frameworks is the so-called hybrid modelling [2,28,34,41,47,48].Hybrid multiscale models are composed by different sub-models for different levels of biological organisation (intracellular processes, cell-to-cell interaction, secretion and transport of signalling cues, etc.), and each of these is modelled in terms of different mathematical descriptions (ODEs, cellular automata, PDEs, etc.).Hybrid models are usually characterised by the use of individual-based models for the dynamics of at least one of the cellular compartments considered in the model [42,43].Other phases (e.g.cellular populations not modelled as individuals and fluid phases such as blood or interstitial fluid) are modelled by means of PDEs as continuous phases [28,34].The individual-based model are often supplemented with models for the behaviour of individual cells in response to cues such as signalling molecules or nutrient depravation [2].The concentration of nutrients or signalling molecules is normally modelled as a continuum field by means of PDEs of the reaction-diffusion type.Hybrid models have been proposed to study different aspects of tumour growth such as response to therapy [2,48], tumour-induced angiogenesis [34,42,43] and evolutionary dynamics of tumour growth [47].
Another possible approach to multiscale modelling is to use multi-phase models [10,33,45,57].In these models, each cellular type is modelled as a different phase.Multi-phase fluid models have been used to analyse different aspects of tumour growth, where each cell type corresponds to a different fluid [10,33,45].Phase-fields models have been recently used to model tumour-induced angiogenesis [57].
In spite of the considerable effort done in the field of multiscale models to tumour growth, there are many aspects of the dynamics of biological tissues which are still poorly explored within this modelling framework.One of them is the effects of noise.Random effects have been included in several hybrid and multiscale models.For example the models of tumour-induced angiogenesis by McDougall et al. [37,38,50], based on earlier work on a hybrid continuum-discrete model by Anderson & Chaplain [5], or the multiscale models of angiogenesis formulated in [42,43] possess an stochastic element, namely, vessel formation is accounted for in terms of a biased random walk model for the movement of tip endothelial cells.However, a general methodology that specifically allows us to incorporate and analyse the effects of noise at the different scales is lacking.As a first step to fill this gap, we propose in this paper a methodology to formulate stochastic multiscale models of cell population dynamics as well as developing asymptotic and numerical methods for their analysis.
In this paper, we aim at formulating stochastic multiscale models of the dynamics of cellular population which account by fluctuations both at the level of intracellular signalling pathways, due to low protein numbers, and at the level cell-population, due to finite population size effects.Hereafter, we will refer to these two sources of noise as molecular noise and cellular noise, respectively.The aim of this paper is to address the issue of noise in multiscale systems in a systematic way.To this end we set up a framework that allows us to formulate and analyse stochastic multiscale models.
Concerning the set up of our model, we will make the same basic assumption as the one done in [2], namely, we will divide the consider three layers which we identify with processes characterised by widely diverse times scales (see Fig. 1 for a schematic representation of our model and the characteristic time scales involved).We consider a model where we couple the dynamics of the concentration of available nutrient (e.g.oxygen), determined by its rate of supply and consumption by the cells, an intracellular layer, where we consider a model of how the concentration of oxygen regulates the rate of progression through the cell-cycle [1,7], and, therefore, also the birth rate, and, last, a cellular layer where we consider the stochastic dynamics of the population of cells.The intracellular and cellular layers are coupled by a model for the oxygen-dependent birth rate, formulated in terms of a mean first passage time problem.
The remainder of this paper is organised as follows.In Section 2, we describe the formulation of our model.We discuss our stochastic model of oxygen-regulated cell-cycle progression and address the formulation of a model for the oxygen-and-age-dependent birth rate as a mean first passage time problem associated to the stochastic cell-cycle dynamics.We then proceed to formulate an age-dependent stochastic birth-death process for the dynamics of the cell population.In Section 3, we present a WKB asymptotic method to find approximate solutions of the corresponding age-dependent Master Equation.Section 4 is devoted to the formulation of an age-dependent τ -leap method that allows us to perform simulations of the multiscale stochastic model.Last, in Section 5, we discuss our results, the limitations of the present approach and directions for future research.

General structure of the stochastic multiscale model
Before going through a detailed discussion of the different elements involved in the formulation of the stochastic multi-scale model, we proceed to a detailed description of the general structure of the model which is closely related to that of the model proposed in [2].
The model we present in this article integrates phenomena characterised by different time scales, as schematically shown in Figure 1, which include oxygen delivery to and consumption by the cell population, the dynamics of the population of cells under the restriction of a oxygen supply at a finite rate, and cell division proliferation and apoptosis.Its structure is therefore quite complex and, for this reason, before presenting the sub-models involved in the description of each of these processes, we explain the overall structure of the modelling framework.The oxygendependent birth rate is modelled in terms of a mean first passage time problem, which is then used in the Master Equation which determines the stochastic dynamics of the cell population.Cells consume oxygen, and therefore their dynamics regulates the concentration of oxygen, thus closing the feed-back loop.We also show the corresponding characteristic time scales: tenths of milliseconds for oxygen, minutes or hours for the intracellular processes, and days for cells.
The approach we use is a natural generalisation of the standard continuous-time birth-and-death Markov process and its description via a Master Equation [17].As we will see, the consideration of the multiscale character of the system, i.e. the inclusion of the physiological structure associated to the cellcycle 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.
Regarding the particulars of each sub-model involved, the model of oxygen delivery is a stochastic differential equation where the oxygen is supplied at a constant rate F and consumed by cells (see Fig. 1).The stochastic character of the equation governing the evolution of the concentration of oxygen arises from the fact that the number of cells at a given time is a stochastic variable.
The second sub-model (i.e. the intracellular model) considered in our multiscale framework is an stochastic model of oxygen-regulated cell-cycle progression (see Fig. 1).This sub-model is formulated using the standard techniques of chemical kinetics modelling [20] so that the mean-field limit of the stochastic model corresponds to the deterministic cell-cycle model formulated in [1].This model provides the cell-cycle status, i.e. the number of molecules of each protein involved in the model from which we derive whether the G1/S transition has occured, for cell of a given age, a.The cell-cycle status of a cell of age a is determined in terms of whether the abundance of certain proteins which activate the cellcycle (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 thresold 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-passege time problem (MFTP) in which one analyses the probability of a Markov process to hit a certain boundary [17].The rate at which our cell-cycle model hits the cyclin activation threshold, i.e. the rate at which cells go through the cell-cycle restriction point, is taken as proportional to the birth rate.The birth rate is a function of the age of the cell as well as the concentration of oxygen, as oxygen abundance regulates the rate of progress through the cell-cycle.
The third and last sub-model (i.e. the cellular model), corresponds to the dynamics of the cell population and it is governed the Master Equation for the probability density function of the number of cells [17].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 age-dependent and 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 Multiscale Master Equation does not present the standard form for unstructured populations, rather, it is an age-dependent Master Equation.
The detailed description involved in each one of the sub-models summarised here is the object of the remainder of Section 2.

Background on cell-cycle modelling
The cell cycle is the sequence of events by which a growing cell duplicates all its components and divides into two daughter cells, each with sufficient machinery and information to repeat the process [4].The cell-cycle is usually divided into four phases: G 1 ; S, G 2 ; and mitosis M .In G 1 (G =gap), the cell is not committed to division and the chromosomes do not replicate.Replication of nuclear DNA occurs during the S phase, whereas completion of mitosis occurs in the final M phase.The interval between DNA replication and division is called the G 2 phase.The gap phases G 1 and G 2 give the cell additional time for growth.The cell also passes through two irreversible transitions.The first of these transitions occurs at the end of G 1 and is called "Start".During G 1 the cell monitors its environment and size.When the external conditions and the size of the cell are suitable, the cell commits itself to DNA synthesis and division.This transition is irreversible: once the cell enters the S phase and DNA replication commences, division has to be completed.The second transition,"Finish", occurs when DNA replication is completed.
Once the cell has checked that DNA and chromatide alignment have occurred, the Finish transition is triggered and the cell finally divides into two daughter cells.A fifth state, the so-called G 0 state, is defined to refer to cells that have abandoned normal progression through the cell-cycle and become quiescent.In this state most (although not all) of the cell functions are suspended, most notably, proliferation.
Cell cycle events are controlled by a network of molecular signals, whose central components are cyclindependent protein kinases (Cdks).In the G 1 state, Cdk activity is low, because its obligate cyclin partners are missing, because cyclin mRNA synthesis is inhibited and cyclin protein is rapidly degraded.At Start, cyclin synthesis is induced and cyclin degradation is inhibited, causing a dramatic rise in Cdk activity, which persists throughout S, G 2 and M. High Cdk activity is needed for DNA replication, chromosome condensation, and spindle assembly.At Finish, a group of proteins, making up the anaphase-promoting complex (APC), is activated [63].The APC attaches a "construction label" to specific target proteins, which are subsequently degraded by the cell's proteolysis machinery.The APC consists of a core complex of about a dozen polypeptides plus two auxiliary proteins, Cdc20 and Cdh1, whose apparent roles (when active) are to recognize specific target proteins and present them to the core complex for labelling [61,63].Activation of Cdc20 at Finish is necessary for degradation of cohesins at anaphase, and for activation of Cdh1.Together, Cdc20 and Cdh1 label cyclins for degradation at telophase, allowing the control system to return to G 1 .We must distinguish between these two different auxiliary proteins, because Cdc20 and Cdh1 are controlled differently by cyclin-Cdk, which activates Cdc20 and inhibits Cdh1.
In [58], Tyson & Novak describe a model for the irreversible transitions "Start" and "Finish" which regulate cell-cycle progression.The model we put forward assumes that these transitions occur by means of bifurcations of the regulatory system which lead to the creation and destruction of stable steady states of the molecular regulatory system of the cell division process.
The dynamics of the cell-cycle can be affected by environmental conditions, in particular, by the level of extracellular oxygen: It is well documented that low oxygen concentrations (hypoxia) alter progression through the cell division cycle [18], and the G 1 /S transition, in particular.In reference [1], it was assumed that the response of this transition to hypoxia is mediated by the protein p27, an element of the Cdk network whose production is upregulated under hypoxia [16,18], although recent studies cast some doubts on the role of p27 as the mediator of hypoxic effects on cell-cycle progression [9,23].In our model, we assume p27 mediates hypoxia-induced arrest of the G 1 /S transition by inhibiting cyclin-Cdk complex formation and, thereby, inhibiting DNA synthesis.
A modification of the model by Tyson & Novak [58] is proposed in [1], where the effects of hypoxia through p27 levels on the cell-cycle are considered.p27 inhibits the formation of cyclin-CDK complex.In turn, p27 levels raise in the presence of hypoxia.The set of ordinary differential equations introduced in [1] to model the effect of hypoxia on the Start transition is the following: where x and y are the concentrations of active Cdh1/APC and concentration of cyclin-CDK complexes, respectively, z, concentration of p27, O 2 oxygen concentration u generic activator, η the cell growth rate, m is the mass of the cell, and m * is the mass of an adult cell.The a i (i = 1, 2, 3, 4), c zi , (i = 1, 2), b i (i = 3, 4) are rate constants and the J 3 and J 4 are Michaelis-Menten constants.Tyson & Novak [58] scale their equations so that the total concentration of Cdh1 (active plus inactive) is normalized to 1 and the Michaelis-Menten constants J 3 and J 4 are such that J 3 ≪ 1 and J 4 ≪ 1.Although in the model proposed in [1], the cyclin considered to be involved in the G1/S transition is CycB and its inhibitor is assumed to be APC/Cdh1, recent work suggests that a more accurate depiction of the situation would involve considering CycE and its inhibitor SCF, instead [60].In fact, for mammals, cyclin D is involved the regulation of the slow dynamics of phase G1 (inhibited by p27), whereas cyclin E regulates the fast dynamics (see reference [54] for a detailed description).A fully accurate acount of the regulation of the G1/S transition in mammalian cells should take the presence of these two cyclins into consideration, rather than lumping their effects into a single compound.

Stochastic formulation
We now proceed to formulate a stochastic model for the oxygen-regulated progression through the cellcycle as a Markov process in terms of a Master Equation.The resulting model will be analysed using large system size, WKB asymptotic [3,30,35].This (sub-)model provides the proliferation rate of the cells as a function of the extracellular oxygen.This information will then be used within the cellularscale population model as a parameter, i.e. the oxygen-age-dependent birth rate in the Master Equation describing the dynamics of the cellular phase.The model we propose here is based on the same basic principles [58] as the one formulated in [1].Tyson & Novak [58] proposed a model for the G 1 /S transition in which the central element of the model is the mutual inhibition between the active form of Cdh1/APC, an inhibitor of cell-cycle progression, and CycB-CDK, whose activity is needed in order for the cell-cycle to undergo the aforementioned transition.This mutual inhibition gives rise to a bistable system with two stable steady-states: the so-called G 1 fixed point, where Cdh1 activity is close to its maximum and CycB activity is virtually non-existent and . Schematic representation of the reactions involved in the stochastic of oxygenregulated cell-cycle progression.M is the mass of the cell, X(X 1 ) is the number of active (inactive) Cdh1 molecules, , Y is the number of cyclin-CDK complexes, and Z is the number of p27 molecules.Reactions (a) and (b) correspond to enzyme-catalysed activation and inactivation of Cdh1, respectively.Note that the inactivation reaction is upregulated by CycB (Y ) and modulated by cell size (M ).Reactions (c) and (d) determine the dynamics of the number of active CycB and p27 molecules.
CycB is synthesised at a constant rate and degraded at a rate which depends on both active Cdh1 and p27.p27 is synthesised at a size-depending rate and degraded at an oxygen-depending rate.According to [1], The rate constants are given in Table 2 the S-G 2 -M fixed point where the opposite is true.Furthermore, Tyson & Novak [58] assume that Cdh1 inhibition by CycB is modulated by cell-size: inhibition is initially poor when cells have just divided and had not reach the necessary critical size to enter the S-phase but it gets enhanced as cells grow and approach such critical size.Mathematically, such cell-size regulations induces a saddle-node bifurcation in which the G 1 fixed point is destroyed when cell size (mass) reaches a critical value, which forces the system to increase CycB activity and to enter the S-phase.In [1] a modification of this very simple was proposed, whereby a further inhibitor of cyclin activity, p27, was introduced.The activity of p27 is known to be upregulated by lack of oxygen (hypoxia), which delays the onset of the G 1 /S transition.This model allows us to couple the rate of cell-cycle progression with the abundance of oxygen and, therefore, to analyse the effects of fluctuations in the supply of oxygen on tumour growth [2].
The reactions involved in our stochastic, oxygen-regulated cell-cycle progression model are schematically shown in Fig. 2. In Fig. 2, M stands for the mass of the cell, X(X 1 ) is the number of active (inactive) Cdh1 molecules, E 1 (E 2 ), the number of Cdh1-activating (inactivating) enzymes and C 1 (C 2 ), the number of X 1 E 1 complexes (XE 2 ).Moreover, Y and Z refer to the number of cyclin-CDK complexes and the number of p27 molecules, respectively.
Reactions Fig. 2(a) and (b) correspond to enzyme-catalysed Cdh1 activation and inactivation reactions, respectively.Note that, as per [58], Cdh1 inactivation is upregulated by active CycB (Y ) and it is modulated by cell growth (M ).Reaction Fig. 2(c) accounts for the dynamics of the number of active CycB molecules, Y : CycB is synthesised at a constant rate and degraded at a rate that depends both on the number of active Cdh1 (X), thus closing the negative feed-back loop of mutual inhibition between Cdh1 and CycB, and also on the number of p27 molecules (Z), which implements in the model the role of p27 as an inhibitor of cyclin activity.Last, reaction Fig. 2(d) determines the dynamics of the number of p27 molecules, Z: p27 is synthesised at a cell-size depending rate and degraded at an oxygen-depending rate in such a way that, when oxygen is scarce, degradation of p27 is down-regulated.This effect yields a build-up of p27 which delays cell-cycle progression by increasing inhibition of cyclin activity.The reader is referred to [1] for full details on the biological rationale of this model.
The stochastic model is thus specified in terms of the state vector, − → X (a): where a stands for age which is taken to be the time elapsed since the last cell division.The dynamics of the model is described by the probability density of the system being in state − → X at age a, Ψ ( − → X , a), whose dynamics is given by the Master Equation: where W i ( − → X , a) are the transition rate corresponding to each of the elementary reactions involved in the model shown in Fig. 2, and r i is a vector whose entries correspond to the increase in the number of molecules of each molecular species when reaction i fires up, i.e.P ( The transition rates corresponding to the enzymatic reaction Fig. 2 are given in Table 1.We have used the law of mass action (LMA) to model the kinetics of the chemical reactions shown in Fig. 2 [20], including the enzymatic reactions Fig. 2(a) and (b).We have chosen LMA kinetics to model these reactions instead of Michaelis-Menten kinetics for technical reasons that have to do the asymptotic analysis of Eq. (2.1) that we carry out in the next section: WKB asymptotics demand the transition rates W i satisfy certain scaling laws (see Eq. (2.3) below).Such scaling relation is not satisfied by the Michaelis-Menten rates, so that we need to resort to LMA kinetics.However, this means that one needs to be careful when parametrising the model, as the parameter values given in [58], where Michaelis-Menten kinetics were used to model enzyme-catalysed activation and inactivation of Cdh1, are not directly applicable to our model.In Appendix A we address this issue and find the relation between the parameters used in [58] and ours.

Model Analysis: WKB approximation
The methodology we use to analyse our model is based on WKB asymptotics and was first proposed by Kubo et al. in [30].In [30], Kubo et al. have proved that, under the appropriate scaling assumption, the time-dependent solution of the Master Equation (ME), Eq. (2.1), can be approximated by a function of the same form as its equilibrium solution, namely the exponential of a homogeneous function, which we call S, of − → X , where Ω is some measure of system size [6].Kubo et al. [30] have shown that the transition rates W ( − → X , r, a) must be homogeneous functions of X to obtain a solution of the ME of the form of equation (2.2), Accordingly, the probability of a given reaction to occur within an infinitesimal interval of time is proportional to the size of the system, Ω, and is determined only by the state of the system, represented by the set of intensive variables − → x .The definition together with equation (2.3), enables us to write the ME (2.1) in WKB from where we have used that e −r•(∂/∂x) is the generator of the translations in the space of states of the system.
To proceed further, we consider the characteristic function of ψ( − → x , a), and its associated cumulant generating function q(u, a) ≡ log(Q(u, a)) [3,30], as the cumulants q n (a) of ψ( − → x , a) can be obtained from the expansion: where u n stands for the n-adic product defined as (u n ) j1,j2,...,jn ≡ n i=1 u ji and "•" denotes full contraction over all of the n indexes.It can be shown [3] that, in terms of the characteristic function Q(u, t), Eq. (2.4) is given by: where w i (v, a) is the Fourier transform of w i ( − → x , a).Kubo et al. [30] showed that Eq. (2.5) is the starting point for an asymptotic expansion of the WKB type, where a closed hierarchy of ordinary differential equations for the cumulants (q n (a)) of the process is obtained.It has been proved by Kubo et al. [30] that, for arbitrary n, the cumulants of the probability distribution (e.g., (q 1 ) i = x i ; (q 2 ) ij = x i x j − x i x j ...) satisfy the following scaling relation: q n (a) = ε n−1 q n1 (a)+ ε n q n2 (a)+O(ε n+1 ), where ε = Ω −1 .This scaling, in turn, yields a consistent asymptotic expansion leading to a system of ordinary differential equations for the cumulants q n (a) in terms of all the cumulants of lower order, q 1 (a), . . ., q n−1 (a).
The above scaling for the order-n cumulants imply that a Gaussian approximation of the process can be obtained such that X(t) = N N (q 11 , N −1/2 q 21 ) where q 11 (a) is the lowest order approximation for the first cumulant (i.e. the first moment, q 1 ), which satisfies the mean-field equations [3,30]: and q 21 (a) is the lowest order approximation for the second cumulant (i.e. the covariance matrix, q 2 (a)), whose components satisfy the set of ODEs: where For a full details about the WKB method, including detailed derivations of Eqs.
(2.6), the reader is referred to [3,30].We use equations (2.6) and (2.7) to formulate the systems of ODEs for the leading-order contributions to the first and second cumulants (i.e. the first and second moments, respectively) By substituting the corresponding values of w( − → x , r, a) and r from Table 1, into Eq.(2.6), where q 11 = − → x = x is mean vector.We obtain the following equation for each element of the mean vector, where abusing of notation ) Similarly, using equation (2.7) we obtain the corresponding set of ODEs which, coupled with Eqs.2.8, allows us to obtain the entries of the covariance matrix, σ(a) = (Q ij (a)).
The Fig. 3 shows the comparison between the numerical solution of Eq. (2.8) and direct numerical simulation of the stochastic system using Gillespie algorithm [20] Finally, using the previous results regarding WKB asymptotics of the Master Equation, we obtain a Gaussian approximation to the solution of Eq. (2.1) is: where d is the number of chemical species in our cell cycle model, |•| represents determinant, () −1 inverse, X(a) is q 11 and σ(a) is the (symmetric) covariance matrix with Q ij are the components.

Modelling the age-dependent birth rate
Activation of many regulatory pathways, in particular, the cell-cycle, depends upon a particular component of the regulatory system reaching a critical activation value.In the case of our stochastic model for cell-cycle progression, cells go through the G 1 /S transition when the level of activity of CycB reaches a threshold value.We assume that, after the cell goes through this transition, the cell completes the cell-cycle and eventually divides after an average time τ p has elapsed.Therefore, given some external conditions (in our case, such conditions are characterised by the concentration of oxygen), the probability of a cell to divide after certain age is equal to the probability that the corresponding protein has reached its critical value.This can be formalised in terms of a first passage time problem [46] whose solution provides precisely this probability and its derivative with respect to the age gives us the birth rate at age a. Once we have obtained the corresponding birth rate in this manner we can use it to parametrise a Master Equation for the stochastic evolution of the cell population.
To calculate the birth rate (i.e. probability of birth per cell and per unit time) in terms of a first passage time problem, we consider the G 1 /S transition occurs when cyclin activation reaches a threshold value and, therefore, we can define the birth rate as a function of the time as b(a We are interested in a region where the y is greater than a constant, k, with all other variables varying over their whole ranges.So, in this case, the expression of G reduces to:  G(a) can be re-written by using the error function as follows: where φ(y) = 1 2 (y − y(a)) 2 σ y (a) −1 .Using the divergence theorem and Taylor's expansion, the first integral in Eq. (2.12) is:

G(a)
whereas the second integral in Eq. (2.12) can be approximated by extending the integration domain to ±∞: We have checked by means of numerical integration that the error introduced by replacing the (y < k)−integration domain by R in I 2 is negligible.
Finally, we find that the ODE G(a) must satisfy is: Fig. 4 shows the probability of cyclin-Cdk exceeds its threshold, 1 − G(a), using Eq.(2.11) coupled with Eq. (2.8) and (2.7), and compares the result with the Gillespie stochastic simulation algorithm [20].
In Fig. 5, we have plotted the birth rate for different oxygen concentrations.As we can see, the peak in b(a) occurs later times for lower oxygen concentrations, this confirms that fact of oxygen delays the occurrence of the G 1 /S transition.And in Fig. 6 we observed the probability of birth has more variance if we increase the noise in the process.

Formulation of the age-structured Master Equation
We are now in position of formulating our model for population dynamics.Using the above framework calculate the birth rate as a function of cell age.We formulate an age-dependent birth-death process.with Ω = 10 7 , O2 = 1.0 and initial condition x = (5, 1.4, 0.01, 0.01, 0, 0.99, 0.01, 0, 0.01).
Parameters such as the oxygen-age-dependent birth rate are determined in term of the models analysed within the intra-cellular scale described in the previous section.To derive the age-dependent Master Equation, we consider the following identity: where W (n(a), a, t) = (ν + b(a))n(a), the birth rate is given by b(a) = − dG(a) da and ν is the constant death rate.
The p 0j births at age a j is provided by a Poisson random variable with mean (and variance) b(a j )n(a j )δt, P(b(a j )n(a j )δt) = p 0j .The probability of p 0 births is the follow: where P Pa j (p 0j ) = e −b(aj )n(aj )δt (b(aj )n(aj )δt) p 0 j p0 j !denote the probability that P(b(a j )n(a j )δt) = p 0j .When a birth even occurs, the number of cell with a = 0 is increased by 2(p 0 ), whereas the number of cells with age a = a j is reduced by p 0j .By re-arranging, Eq. (3.1) take the limit δt → 0, we obtain: ∂P(n,a,t) ∂t + ∂P(n,a,t) ∂a = W (n(a) + 1, a, t)P(n(a) + 1, a, t) −W (n(a), a, t)P(n(a), a, t). (3.3) In order to find an approximate a solution of Eq. (3.3) we apply the WKB method put forward by Kubo et al. [30] have shown that the transition rates W (n, a, t) must be homogeneous functions of n to obtain a solution of the ME of the form P (n(a), a, t) = C exp(−N s(n(a), a, t)) with n(a) = n(a)/N W (n(a), a, t) = N w(n(a), a, t). (3.4) Accordingly, the probability of a given reaction to occur within an infinitesimal interval of time is proportional to the size of the system and is determined only by the state of the system and is determined only by the state of the system, represented by the set of intensive variable n(a).The definition P (n(a), a, t) = N P(n(a), a, t), together with equation (3.4), the Eq. ( 3.3) has the following expression: where we have used that e −(∂/∂n(a)) is the generator of the translations in the space of states of the system.
Let us consider the cumulant-generating function defined as q(u, a, t) = ln (Q(u, a, t)), where Q(u, a, t) is the characteristic function of P (n, a, t), defined as its Fourier transform: The cumulants, q k (a, t), can be obtained from q(u, a, t) as the coefficients of the expansion: The procedure we follow next, which is based on the work of Kubo et al [30], is to write down an equation for Q(u, a, t) and construct asymptotic approximations for Q(u, a, t), q(u, a, t) and q k (a, t).These approximations will provide systems of ODEs for the cumulants and for the moments of the solution of the ME.We are interested to the system of equation for the first moments (q 1 (a, t)) and for the elements of the covariance matrix (q 2 (a, t)), which will enable us to study the average behaviour of the system and the Gaussian fluctuations around it.
To obtain the equation for Q(u, a, t), we take the Fourier transform of the ME Eq. (3.5): 1 N ∂Q(u,a,t) ∂t The integral on the right hand side of Eq. (3.7) e i u n w(n, a, t)P (n, a, t) dn. (3.8) Furthermore, substituting Eq. (3.8) in Eq. (3.7), rearranging terms, and recalling that the Fourier transform of the product of two functions equals the convolution of the corresponding Fourier transforms, we finally obtain: where ω(v, a, t) is a Fourier transform of w(n, a, t).
Eq. (3.9) needs to be supplemented with boundary conditions for P (n, a = 0, t), i.e. for the probability of the number of births within the time interval (t, t + δt) to be p 0 = n/2.The number of births within this time interval can be expressed as the sum of all the births at time t within each age-group a j , p 0j , i.e. p 0 = j∈I(t) p 0j .Therefore, where p 0j = Y j (b(a j )n(a j )δt) is a Poisson-distributed random variable with parameter λ j = b(a j )n(a j )δt, i.e.
and, consequently, the corresponding cumulant-generating function for the random variable p 0j is Q j (s) = e λj (e is −1) .Now, since P (p 0 ) is a convolution, the cumulant generating function for P (p 0 ), Q(s, a = 0, t) is given by [24]: which, taking the limit δt → 0, results in:
Let us focus on the first cumulant q 1 (a, t) = n(a, t) , which is given by the O(ǫ) terms in Eq. (3.12).Expanding the exponentials within the right hand side of Eq. (3.12) and keeping only first order terms we obtain: where ǫ ≡ N −1 and m k (ǫ, v, a, t) is defined by: exp Eq. (3.13) needs to be properly balanced with respect to the small parameter ǫ.Such balance is achieved if, to leading order, m 0 (ǫ, v, a, t) = O(ǫ 0 ), which is, in fact, satisfied as m 0 (ǫ, v, a, t) = 1 (see Eq. (3.14)), and q 1 (a, t) = O(ǫ 0 ).
To leading order, the equation for q 1 (a, t) reads: e q(−v,a,t) ω(v, a, t) dv, (3.16) whereas the corresponding equation for q 2 (a, t) is given by: e q(−v,a,t) ω(v, a, t). (3.17) The leading order approximation of Eq. (3.17) is obtained as follows.From the definitions of the quantities h j and m j we obtain the following expression for h j : For j = 1, from Eq. (3.14), we have that m 1 = h 1 .Taking into account the scaling substitution for m k and q k , we obtain the following O(ǫ) approximation: Thus, substituting Eq. (3.18) and keeping only the leading order contributions, the equation for q 21 (a, t) reads: (−ivq 21 (a, t) + 1) e q(−v,a,t) ω(v, a, t).
Furthermore, the scaling substitution for q k (a, t) and Eq.(3.6) lead to: This almost completes our derivation of the evolution equations for q 1 (a, t) and q 2 (a, t) to leading order Eq.(3.19) can be rewritten as ) we obtain that, to leading order, the equations for q 11 (a, t) and q 21 (a, t) are the following equations: where w(v, a, t) is the Fourier transform of w(n, a, t) so the quantity c(q 11 , t) = − (ν + b(a)) q 11 .Finally, the evolution equation for n(a, t) = q 11 (a, t) is given by: ∂n whereas the equation for the variance and σ n (a, t) = q 21 (a, t) is given by: ∂σn(a,t) dt

Gaussian approximation for the stochastic multiscale model
After the derivation of the previous sections, we can now write down a Gaussian approximation to the multiscale model schematically represented in Fig. 1.The coupled system which determines such approximation is as follows: -Oxygen concentration, O 2 : where O 2 is the oxygen concentration, which is consumed by the total population, M (t), S is the rate of delivery of oxygen to the population, which we assume to be constant, and κ is the per-cell oxygen consumption rate, which we assume to be age-independent -Age-and oxygen-dependent birth rate (Section 3): -Mean and variance of the age-structured population (Sections 3.1 and 3.2): )

Numerical scheme: Age-dependent τ −leaping method
In this section we describe a stochastic numerical method to compare the solution of ME (3.3)-(3.2) with the solutions of WKB approximation given by system of equations (3.24-3.26).
We use an extension of the τ -leaping's simulation method in the birth-death process described by ME (3.3)-(3.2),where the probability of birth is described by Eq. (2.11).Using τ -leaping in every age group present in the population, n(a i ), we generate a time step in each age, τ i .The global time step, τ , is chosen to be the minimum of τ i ∀i.At this point, we can calculate the number of of birth and death events that will occur within time span τ at every age group n(a i ), which are given by Poisson deviates P(d ij τ ) where: where Eqs.(4.1) and (4.2) are the probabilities of birth and death during a (short) time interval of duration τ for cells of age a(i).The events whose probabilities are given by Eqs.(4.1) and (4.2) have associated changes in the population n(a(i)) given by r i1 = −1 and r i2 = −1, respectively.In addition, a birth event (whose probability is given by Eq. (4.1)) have associated a change in the new born cells, i.e. cell with age a = 0, given by r 0 = 2.We first introduce the method and how to generalise it to our age-structured problem.We then proceed to describe the algorithm we use in our simulations.

Age-structured τ -leaping method
In this section, we are going to explain the τ -leaping method describing in [12,21] and the modification that we have carried out.Gillespie's stochastic simulation algorithm (SSA) [20] provides exact sample paths of the probability density which solves the corresponding Master Equation.It does so by carrying out every single event in the sample path which means that this method can be painfully slow in particular for systems involving both fast and slow processes.The SSA also generates a vast amount of detailed information on a particular sample path whose knowledge may be unnecessary in many applications.In view of this situation, the τ -leaping method [21] was proposed which aims to speed up the simulation of Markov processes by answering the following question: How often does each process happen in the next specified time interval τ ?The mathematical foundation of this algorithm can be traced back to the so-called representation of the Markov processes.Specifically, assume a Markov process, X(t), whose probability density satisfies a Master Equation: Each reaction, i, involved in the process described by Eq. ( 4.3) is characterised by a exponential waiting time distribution [20]: We therefore know that the number of times channel i fires in the time interval (0, t) is given by P t 0 W i (X, s)ds where P(λ) is a random number sampled from a Poisson distribution with parameter λ [24].Therefore, X(t) can be formally represented as: Using the Markov property, X(t + τ ) can be expressed in terms of X(t) as: Thus, the question that the τ -leaping method aims to address can be reformulated into the following form: Under which conditions Eq. (4.4) can be approximated by: In other words, we need to estimate τ so that t+τ t W i (X, s)ds ≃ W i (X, t)τ .This is a problem for which no systematic solution has been devised so far.Instead, estimates, referred to as leap conditions, have been derived based on different assumptions [12,21].Roughly speaking, τ needs to be chosen so that |W i (X, t) + k i ) − W i (X, t)| is bounded by a small quantity.
In order to generalise this method to our age-structured model, we first consider an age-number vector whose components are n(a i ) with i ∈ I t where I t spans the set of ages actually represented in the population, i.e. such that n(a i ) > 0 at time t.The first issue is to choose a value for τ that satisfies a leap condition.The simplest generalisation of the standard τ -leaping algorithm is to choose a leap value, τ i , for every age a i with i ∈ I t , according to the criteria established in [12] and then choose the minimum among all of them.We start by computing the time τ i for every age, a i in i ∈ I t , using the method formulated by Cao et al. [12], which improves upon previous results by Gillespie & Petzold [22], as it adheres more closely to the leap condition by uniformly bounding the relative changes in the propensity functions.Moreover, this method is faster than the procedure suggested in [22], since the number of auxiliary steps used to calculate the time leap increases linearly with the number of reactant species, rather than quadratically with the number of reaction channels [12].
The underlying strategy of this τ i -selection procedure is to bound the relative changes in the number of individuals of populations in such a way that the relative changes in the propensity functions will all be approximately bounded by a specified value ǫ (0 < ǫ < 1).Let The τ selection we shall base it on the condition where J rs denotes the set of indices of all population with different age (so j ∈ J rs if and only if d ij is an argument of at least one propensity function).
Since the Poisson random variables P i (d ij τ i ), on the right-hand side of Eq. (4.5) are statistically independent and have means and variances d ij τ i , the mean and variance of the linear combination Eq. (4.5) can be computed straightforwardly, Using the same reasoning that was used in deriving the Gillespie and Petzold τ -selection procedure [22], we may consider the bound (4.6) on ∆ τi n(a(i)) to be substantially satisfied if it is simultaneously satisfied by the absolute mean and the standard deviation of ∆ τi n(a(i)): Substituting formula (4.7) into conditions (4.8), we obtain the following bounds on τ so taking At this point, we have computed a time leap for every age, τ i .To take the system forward in time, we take the minimum τ i which satisfies the corresponding leap condition τ = min i∈I {τ i } and calculate the quantities k ij = P(d ij , τ ) and λ i = j k ij .Finally, we make effective the leap by replacing t by t + τ and n(a(i)) by n(a(i)) − λ i .If k i1 = P(d i1 , τ ) is not zero for some age i, it means that in the population with age a i there have been births, so the zero-age population, n(0, t + τ ), is given by n(0, t + τ ) = 2 i P(d i1 , τ ).This also implies the age-number vector has now a new component, i.e.
To calculate d i1 , we need to know the value of the birth rate, b(a i ) at new time t + τ .For this, we resolve the system (3.24-3.25) in the interval of time (t, t + τ ) for each age using the vector − −− → n(a i ), we use implicit 2nd order Runge-Kutta implemented in the GNU Scientific Library (GSL).Now, we have the probability of birth at time t using Eq.(2.11), we obtain 1 − G(a i ) t at time t.We calculate the rate birth for each age,

Description of the algorithm
The simulation algorithm is: 1. Initialisation.Start with an initial age-number vector, − −−− → n(a, 0) at t = 0, whose components, n(a i , 0), are the initial population corresponding to each age group present in the initial population.The dimension of this vector is denoted by I, i.e. the number of age groups with non-zero initial population.At t = 0, b(a i , t = 0) = 0 for all age.2. Compute birth and death rates for each age: d i1 = b(a i )n(a(i)) and d i2 = νn(a(i)).4. Update the system variables: 5. Determine the number of births: n(0, t) = 2 i P(d i1 τ ).If the number of births n(0, t) = 0, we must add the new component to the age-distribution vector, i.e.I ← I + 1. 6. Solve equation for the oxygen concentration in the interval (t, t, +τ ): is the total population.
7. Calculate the birth rate, b(a i ), using the probability of birth at time t using Eq.(2.11), we obtain 1 − G(a i ) t at time t.We calculate the rate birth in each age, b(a i ) = (1−G(ai)t−τ )−(1−G(ai)t) τ .8. Loop through items 2 to 7 until t ≥ T .Simulation results using this algorithm are shown in Figs.7 and 9, where the age-structured τleaping simulation algorithm is compared with the numerical solution of the system of equations Eqs.(3.22) in order to validate the WKB-Gaussian approximation.The Gaussian approximation Eqs.(3.22) implies that a 68.2% of the realisations should stay within the region delimited by the upper and lower boundaries given by n(t, a) ± Ω −α0/2 σ n (t, a). Figure 7 shows that the stochastic sample paths generated by our stochastic algorithm lay between the bounds predicted by the Gaussian approximation.

Conclusions & discussion
In this paper, we present a formulation of a stochastic multiscale model of cell population dynamics where different levels of biological organisation, characterised by different time scales (see Fig. 1) are coupled.We describe the stochastic cellular population dynamics (cellular scale, as per the terminology of [2]) by means of a birth and death process where the birth rate is determined by a model of oxygen-regulated cell cycle progression (intracellular scale, as per the terminology of [2]), thus coupling the intracellular and cellular scales.The resulting coupled system is a stochastic age-structure cell population.The birth rate is described in terms a novel method on a mean-first passage time approach, which allows us to determine the birth rate as the average rate at which key proteins involved in cell-cycle regulation cross a threshold of activation as a function of cell age (i.e.time elapsed since last cell division).Therefore our methodology allows us to analyse the coupling between fluctuations due to molecular noise at the intracellular level with cellular noise due to finite population size.In order to analyse the behaviour of our model, in Section 3, we have extended a WKB asymptotic method initially proposed by Kubo and co-workers [3,30] to obtain an approximate solution of our age structure birth-and-death Master Equation, Eq. (3.5).Carrying out this approximation to lowest order, we obtain a Gaussian approximation to the (age-structured) probability density where the corresponding mean, n(a, t), and variance σ n (a, t) are the solutions of two coupled semi-linear PDEs (see Eqs. (3.24-3.26)).
The semi-linear PDEs for n(a, t) and σ n (a, t) are coupled to the intracellular scale through its dependence upon the age-dependent birth rate b(a, t).This quantity is obtained in terms of a mean-first passage time for the oxygen-dependent cell-cycle model (details are given in Sec.2): We consider that the cell enters the proliferating phase of the cell-cycle when the concentration of CycB exceeds a given threshold.The birth rate is therefore given by the rate at which the system crosses said threshold.In order to solve this problem we carry out a WKB approximation [3,30] to the corresponding Master Equation (se Eq. (2.1) and Table 1) to obtain a Gaussian approximation, which is determined by the solution of the system of ODEs for the mean, Eq. (2.6) and for the variance (2.7).This, in turn, allows us to find an approximation to the solution of the mean-first passage time problem and, therefore, to the age-dependent birth rate.
Our formulation for the oxygen-regulated cell-cycle progression is described by the mean-first passage time activation.The birth approach is calculated by the probability of cyclin-Cdk is down a threshold, it is the probability of extinction, Eq. (2.10), so the birth process happens with the complementary probability.To calculate this activation process, we have used WKB method to obtain an approximate solution to the cell-cycle ME.This result is a series of systems of ODEs for the cumulants of order n, we have used order 2 to obtain a solution of the ME, Eq. (2.9).According to the results from our reduced model compare with the Gillespie simulations, the approach is good approximation only for less intracellular noise, big size of the system form order Ω = 10 7 , see Fig. 3, to obtain a better approximation we must use a greater order than 2 in the WKB approximation.The system is coupled through a negative feed-back loop where the population consumes oxygen which, in turn, regulates cell-cycle progression and therefore birth rate.
In order to validate our asymptotic approach in the stochastic multiscale model, we have introduced an age-structured extension of the τ -leaping simulation method in birth-death process associated to Eq. (3.3)-(3.2) with the birth rate calculated by Eq. 2.11.We compare the results of stochastic simulations with our WKB approximation of the age-structured cell population dynamics in Figs. 7, where we show realisations of the process generated by using our age-structured τ -leaping model with the Gaussian asymptotic approximation, 8, where we plot the time evolution of the oxygen concentration, and 9, where we show the number of births for the asymptotic approximation and the stochastic simulations.This comparison allows us to assess both the accuracy and range of applicability of our asymptotic method.Our results confirm that, for the intracellular model considered in this paper, the asymptotic approximation is accurate only if we consider the scale for protein numbers in the intracellular cell-cycle model, Ω, to be large enough.In this case, as shown in Figs. 7, 8, and 9 we find that the agreement between numeric and asymptotic results is good.This limitation regarding the number of proteins stems from the fact that, in the cell-cycle model considered in this paper, the transition which heralds the onset of proliferation (i.e.CycB reaching a threshold value) is shortly preceded from a saddle-node bifurcation.Close to such dynamical phase transition, our Gaussian approximation is expected to be unable to capture the statistics of the fluctuations [8].Thus in order to capture the dynamics of the system with accuracy we have to move well into the small fluctuations regime (i.e.large Ω).This means that our model should be able to capture the behaviour of systems exhibiting no such bifurcations in a much larger regime.
An alternative to the WKB method used in this manuscript would be the use of Van Kampen's systemsize expansion [59], also known as the linear noise approximation (LNA), which is closely related to Kubo et al. [30] WKB approximation.In spite of the similarities between both methodologies, they are not exactly equivalent.The LNA assumes that the random variable can be written as X(t) = φ(t) + Ω 1/2 ξ(t), where φ(t) is the solution of the mean-field equations and ξ(t) is a random variables which satisfies a linear Fokker-Planck equations from, where (linear) ordinary differential equations for the first and second moments of ξ(t) can be derived.The WKB approximation is based on an Ansatz regarding the form of the solution of the Master Equation, which means that, as far as our results are concerned, the equations for the variance of X(t) are not, in general, the same, although the same mean-field approximation is obtained in both cases.
The methodology presented in this paper can be generalised in a number of ways as well as being applied to a variety of situations.One such situation that would be interesting to analyse, in particular in relation to cancer evolution and therapy, is to introduce hypoxia-induced quiescence.Cancer cells have the ability to go into quiescence under hypoxic conditions [1,7], so it would be interesting to factor this capability into our model, specially in the context of competition models between cells which do not undergo hypoxia-induced quiescence (normal cells) and cells which enter quiescence (cancer cells) [2], in particular in the context of predicting the outcome of cell-cycle-specific therapies, where cells are killed upon entering specific cell-cycle stages [44].This analysis is postponed for future work.
A related issue that would be interesting to explore using the methodology put forward in this manuscript is the issue of cancer cell metabolism (i.e. the Warburg effect).One of the several, more remarkable phenotypic changes induced by the Warburg effect, together the an increase of the extracellular acidification rate, is the reduction of the oxygen consumption rate (OCR) [13].It would be interesting to analyse the effect of considering different cell types characterised by different metabolic states, i.e. different OCRs, in particular in its relation to the effects of hypoxia.This issue will be further explored in future research.
Further future work, we propose to extend our stochastic multiscale framework to include spatial degrees of freedom and the coupling to a model of tumour induced angiogenesis.To do this, we need to resort and adapt more efficient methods of stochastic simulations such as the Next Reaction Method [19].

A. Appendix A
This Appendix is devoted to establish the relation between the parameters of our stochastic model, d 1 , d −1 , d 2 d 3 , d −3 , d 4 , and the values of the macroscopic parameters given in [1].In order to achieve this, we apply the quasi-steady state approximation (QSSA) to the reaction scheme: This mean-field equations for this reaction mechanism is given by the system of ODEs Eq. (A.1) for the rates of change of each species involved. .Moreover, we use the conservation laws: dei da + dci da = 0, whereby the total enzyme concentration is constant c i + e i = e 0 (i = 1, 2), and dx1 da + dc1 da + dx da + dc2 da = 0, so that the total concentration of Cdh1 is conserved: In order to proceed further, we introduce rescaled variables:  3,4).This equivalence allows us to give values to the parameters of our stochastic model so that the mean-field behaviour is the same as in the models by Alarcón et al. [1] and Tyson & Novak [58].

Figure 1 .
Figure1.Schematic representation of our multiscale model: Oxygen modulates the progression of cells through the cell-cycle, or, equivalently, their birth rate.The oxygendependent birth rate is modelled in terms of a mean first passage time problem, which is then used in the Master Equation which determines the stochastic dynamics of the cell population.Cells consume oxygen, and therefore their dynamics regulates the concentration of oxygen, thus closing the feed-back loop.We also show the corresponding characteristic time scales: tenths of milliseconds for oxygen, minutes or hours for the intracellular processes, and days for cells.
is the probability of cyclin-Cdk exceeds its threshold value at age a. G(a) = R Ψ ( − → X (a)) d − →X is the so-called survival probability, i. e. the probability that − → X (a) ∈ R. If Ψ is approximated by Eq. (2.9) we can obtain a closed equation for G in terms of G, x and Q ij .
X 1 inactivation of Cdh1/APC where the parameters d 1 , d −1 , d 2 d 3 , d −3 , d 4 are related to those in [1] through the QSSA proposed by Briggs and Haldane.

Table 2 .
Parameters values