A Multi-stage Representation of Cell Proliferation as a Markov Process

The stochastic simulation algorithm commonly known as Gillespie’s algorithm (originally derived for modelling well-mixed systems of chemical reactions) is now used ubiquitously in the modelling of biological processes in which stochastic effects play an important role. In well-mixed scenarios at the sub-cellular level it is often reasonable to assume that times between successive reaction/interaction events are exponentially distributed and can be appropriately modelled as a Markov process and hence simulated by the Gillespie algorithm. However, Gillespie’s algorithm is routinely applied to model biological systems for which it was never intended. In particular, processes in which cell proliferation is important (e.g. embryonic development, cancer formation) should not be simulated naively using the Gillespie algorithm since the history-dependent nature of the cell cycle breaks the Markov process. The variance in experimentally measured cell cycle times is far less than in an exponential cell cycle time distribution with the same mean. Here we suggest a method of modelling the cell cycle that restores the memoryless property to the system and is therefore consistent with simulation via the Gillespie algorithm. By breaking the cell cycle into a number of independent exponentially distributed stages, we can restore the Markov property at the same time as more accurately approximating the appropriate cell cycle time distributions. The consequences of our revised mathematical model are explored analytically as far as possible. We demonstrate the importance of employing the correct cell cycle time distribution by recapitulating the results from two models incorporating cellular proliferation (one spatial and one non-spatial) and demonstrating that changing the cell cycle time distribution makes quantitative and qualitative differences to the outcome of the models. Our adaptation will allow modellers and experimentalists alike to appropriately represent cellular proliferation—vital to the accurate modelling of many biological processes—whilst still being able to take advantage of the power and efficiency of the popular Gillespie algorithm.


Introduction
In a well-mixed solution of chemicals undergoing reactions, non-reactive collisions occur far more often than reactive collisions allowing us to neglect the fast dynamics of motion.We can thus assume that the time between reactive collision events is exponentially distributed with rates which are a combinatorial function of the numbers of available reactants (Gillespie, 1976(Gillespie, , 1977)).This premise is the basis of the Gillespie (1976) stochastic simulation algorithm 1 .The Gillespie algorithm has become a ubiquitous algorithm for the simulation of stochastic systems in the biological sciences, in particular in computational systems biology (Szekely and Burrage, 2014).
However, the Gillespie algorithm is often used inappropriately to represent processes for which the inter-event time is not exponentially distributed.One prevalent example of this is in the simulation of the cell cycle (Baar et al., 2016;Castellanos-Moreno et al., 2014;Figueredo et al., 1 Although Gillespie was amongst the first to popularise the stochastic simulation algorithm and derived it from physical considerations, its conception dates back to the work of Doob (1945) and Feller (1940).It was derived in a different form by Kurtz (1972).
A multi-stage representation of cell proliferation as a Markov process 3 2014; Mort et al., 2016;Ryser et al., 2016;Turner et al., 2009;Zaider and Minerbo, 2000) (see figure 1).The assumption of memorylessness, and consequently exponentially-distributed cell cycle times, means that with high probability a daughter cell may divide immediately after the division event which created it.This is not biologically plausible since each cell is required to pass through the G 1 , S, G 2 and M phases of the cell cycle before division, and these phases (in particular S-phase) are rate-limiting.Probability density (× 10 -3 ) data exponential Fig. 1 The poor agreement between the best-fit exponential distribution (red curve) and experimentally determined cell cycle times in NIH 3T3 mouse embryonic fibroblasts grown in vitro (grey histograms).
The rate parameter of the exponential distribution was fitted by minimising the sum of squared residuals between the curve and the bars of the histogram.
As an illustrative example, Turner et al. (2009) present a well-mixed stochastic model of small populations of cancer stem cells, which they use to suggest treatment strategies.Both analytical and simulated results pertaining to the survival of the tumour are based on the assumption that inter-division times are exponentially distributed.We demonstrate later in this paper that using the correct cell cycle time distributions (CCTDs) could alter their results leading to different suggested treatment strategies.
Similarly, in a spatially extended context, Baker and Simpson (2010) investigate the role of spatial correlations on individual-based models of cell migration, proliferation and death, designed to represent experimental assays of cell behaviour in culture.The individual-based models they employ are on-lattice, volume-exclusion processes in two and three dimensions in which cell migration, proliferation and death events are all considered to be exponentially distributed.Cell proliferation occurs when a cell is chosen to divide, placing a daughter cell at one of its neighbouring lattice sites.Baker and Simpson (2010) demonstrate the effects of spatial correlations in these models.In particular, they find that changing the motility of the cells can alter their net rate of growth in comparison to the logistic growth predicted by a simple mean field assumption which neglects the effects of correlations.This effect is due to the fact that cell motility serves to break up correlations allowing more proliferation events to occur in comparison to the lower motility case.By explicitly considering the correlations between the occupancies of pairs of lattice sites, Baker and Simpson (2010) derive a more accurate population-level model which better represents the growth in the number of cells over time for a diverse range of parameter values.We will investigate the effects of incorporating more realistic CCTDs on the outcomes of the model simulations.
Non-Markovian simulation methods exist for events which do not have exponentially distributed inter-event times (Boguná et al., 2014).However, these algorithms are often difficult to understand and complex to encode since we are required to keep track of every cell individually.This presents a potential barrier to their use and consequently a barrier to the appropriate modelling of CCTDs.
Given the ubiquity of the Gillespie algorithm, it would be significantly more beneficial if we could decompose the cell cycle into a series of exponentially distributed events which could be naturally encoded in the framework of the Gillespie algorithm.One potential solution to this problem is the use of the hypoexponential family of distributions.It has been suggested that these distributions can be used to accurately represent phases of the cell cycle (and, by closedness of the sums of these distributions, the cell cycle itself) (Stewart, 2009).Hypoexponential distributions are made up of a series of k independent exponential distributions, each with its own rate, λ i , in series.If k is large then these models may face issues of parameter identifiability.
Recently, Weber et al. (2014) have suggested that a delayed hypoexponential distribution (consisting of three delayed exponential distributions in series) could be used to appropriately represent the cell cycle.These delayed exponential distributions represents the G 1 , S and a combined G 2 ßM phases of the cell cycle.Their model is an extension of the seminal stochastic cell cycle model of Smith and Martin (1973) who use a single delayed exponential distribution to capture the vari-ance in the cell cycle.Delayed hypoexponential distributions representing periods of the cell cycle have been justified by appealing to the work of Bel et al. (2009).Bel et al. (2009) showed that the completion time for a large class of complex theoretical biochemical systems, including DNA synthesis and repair, protein translation and molecular transport, can be well approximated by either deterministic or exponential distributions.
In this paper we consider two special cases of the general hypoexponential distribution: the Erlang and exponentially modified Erlang distribution which, in turn, are special cases of the Gamma and exponentially-modified Gamma distributions.For reference their PDFs P E and P EME , respectively, are given below: (1) With suitable parameter choices, both distributions have been shown to provide good fits to large numbers of experimentally-derived CCTDs (Golubev, 2016)  As special cases of the hypoexponential distributions, these distributions also have the significant advantage that they can be simulated using the ubiquitous Gillespie stochastic simulation algorithm.This will allow for the appropriate representation of CCTDs in stochastic models of cell populations, in contrast to the inappropriate exponentially distributed times which are commonly used (Baar et al., 2016;Castellanos-Moreno et al., 2014;Figueredo et al., 2014;Mort et al., 2016;Ryser et al., 2016;Turner et al., 2009;Zaider and Minerbo, 2000).Additionally the two and three parameters (respectively) of the Erlang and exponentially modified Erlang distributions (respectively) simplify parameter identification in comparison to more highly parametrised distributions.
These two choices (Erlang and exponentially modified Erlang distributions) are not the only nonmonotone distributions which could be used to appropriately represent the cell cycle.However, they are the general, non-monotone, hypoexponential distributions with the fewest number of parameters (two for Erlang and three for exponentially modified Erlang).These features will aid parameter identifiability (few parameters) and crucially mean the distributions can be simulated using the Gillespie algorithm (hypoexponentiality), making these the most suitable distributions to consider.Erlang distributions with the experimental data in comparison to the exponential distribution (c.f.figure 1).In each case the parameters of the distributions are fitted by minimising the sum of squared residuals, Σ, between the curve and the bars of the histogram2 .Clearly, for the exponential distribution (see figure 1), the shape of the curve is incorrect.Consequently the exponential distribution gives a poor representation of the true CCTD with a sum of squared residuals Σ 1.86 ¢ 10 ¡6 .Evidently the Erlang distribution with parameters λ 0.0083 and k 12 gives a much better agreement to the experimental data (see figure 2 (a)), with a minimised sum of squared residuals, Σ 1.23 ¢ 10 ¡7 .Finally, the exponentially modified Erlang distribution with parameters λ 1 0.0251, λ 2 0.0019 and k 26 gives an even better agreement to the data 3 with a minimised sum of squared residuals, Σ 6.01 ¢ 10 ¡8 .The exponentially modified Erlang distribution achieves a minimised sum of squared residuals which is around half that of the Erlang distribution.Never-the-less, both Erlang and exponentially modified Erlang are good candidates for fitting cell cycle time data and can both be simulated within the existing Gillespie framework, so will be considered here.
In Section 2 we begin by outlining a general hypoexponential model of the cell cycle and noting that many previous models of the cell cycle are special cases.By simplifying the model further we demonstrate that the Erlang and exponentially modified Erlang distributions are also special cases.In Section 3, we consider the special case of the Erlang distributed CCTD in more detail.
Undertaking some simple analysis we derive the expected behaviour of the mean cell number in the case of Erlang CCTDs and demonstrate analytically that significant differences can arise in comparison to models in which exponentially distributed CCTDs are used.In Section 4 we demonstrate the utility of our new CCTD representation in stochastic simulations of two biological models in which cellular proliferation is of critical importance.In each case we show, through simulation, that there are important quantitative and qualitative differences between models which represent cell cycle times appropriately and those which do not.We conclude in Section 5 with a short discussion on the implications of our findings.

Multi-stage model of the cell cycle
We divide the cell cycle (with mean length C) into k stages 4 .The time to progress through each of these stages is exponentially distributed with mean µ i .We can represent the progression through these stages of the cell cycle as the following chain of 'reactions' the value of k the more "reaction channels" the Gillespie algorithm must account for which, if rate-limiting, will significantly reduce the algorithm's performance. 4Note, these stages do not necessarily correspond to the traditional (G1, S, G2 and M ) phases of the cell cycle.Indeed, there is good evidence that at least one of the phases is not exponentially distributed (Hahn et al., 2009).Rather they are arbitrary division of the cell cycle which will allow us to recreate the correct CCTD.
where λ i 1ßµ i .
The CCT under this model is hypoexponentially distributed.Although there is no simple closed form for the probability density function of the hypoexponential distribution we can find simple expressions for its mean and variance.The mean is given by the sum of the means of the exponentially distributed stage times k i 1 µ i C and the variance is the sum of the variances of these stage times, By increasing the number of exponentially distributed stages, whilst decreasing their mean duration (in order to maintain the correct mean CCT), we can arbitrarily decrease the variance of the CCT.Many multi-stage models of the cell cycle are special cases of this general model (Golubev, 2016;Hawkins et al., 2007;Hillen et al., 2010;Hoel and Crump, 1974;Kendall, 1948;Leander et al., 2014;León et al., 2004;Nakaoka and Inaba, 2014;Powell, 1955;Smith and Martin, 1973;Weber et al., 2014;Zilman et al., 2010).
We can analyse the cell cycle reaction chain (2) further by considering the associated probability master equation (PME).Let P Ôx 1 , x 2 , . . ., x k , tÕ be shorthand for the probability that there are x 1 cells in stage one, x 2 in stage two and so on.The PME is dP Ôx 1 , x 2 , . . ., x k , tÕ (3) By multiplying the PME by x j and summing over the state space we can find the evolution of the mean number of cells, M j x x j P , in each stage, where x is shorthand for x1 1 ¤ ¤ ¤ x k 1 and P is shorthand for P Ôx 1 , x 2 , . . ., x k , tÕ.Upon simplification we find the following evolution equations for the mean number of cells in each stage (4)

Identical rates of progression
The hypoexponential model's generality is also a significant drawback since it hampers parameter identifiability (Weber et al., 2014).As such, we seek to reduce the number of free parameters in the model whilst maintaining its ability to accurately represent CCTDs.Several authors have suggested using the Gamma distribution to model CCTDs (Hawkins et al., 2007;Kendall, 1948;Nakaoka and Inaba, 2014;Zilman et al., 2010).If we assume that all transition rates, λ i , are identically equal to λ 1 (for i 1, . . ., k) in our general hypoexponential model, then the time to progress through the whole cell cycle is distributed according to the sum of k identically exponentially distributed random variables.It is straightforward to show (using moment generating functions or convolutions) that the CCTD, is Erlang distributed with scale parameter µ Cßk and shape parameter k.In analogy with the general hypoexponential case, if we decrease µ and simultaneously increase k so that µk C remains constant, the Erlang distribution approaches the Dirac delta function centred on C, demonstrating that we can still arbitrarily reduce the variance to match the distribution we are trying to model.

Analysis of the CCDT with equal rates of progression
We now analyse this CCTD model with identical rates of progression, noting that Kendall (1948) studied this case extensively and we draw on some of his analyses below.Although for this special case it is possible to derive a closed form first order partial differential equation for the evolution of the generating function corresponding the the master equation ( 3), solving the associated characteristic equations is analytically intractable for all but the simplest case (k 1) (Kendall, 1948).
Instead we will focus on the mean behaviour of the cell population with which we can make some analytical progress.
In the equal rates case we can sum the individual equations in system (4) to give where Consider the naive one stage (i.e.k 1) cell-cycle mode with mean cell cycle time C: The evolution of the mean number of cells is given by the special case of equation ( 5): dM dt M ßC. (7) In the multi-stage model, under the assumption that all cells are evenly distributed between the stages (i.e.M i M ßk), we can replace M k with M ßk in equation ( 5) to give a closed equation for the evolution of the total number of cells which matches equation ( 7): However, the assumption on the even distributions of cells between stages is incorrect.This leads to differences not just, as might be expected, between the variation exhibited by the multistage and single-stage models, but also between their mean behaviour.In figure 3 (a) a clear difference between the k 1 and k 4 models is evident.The mean total cell number grows significantly more slowly in the k 4 case than the k 1 case.This is true for all models in which k 1. Intuitively, exponentially distributed CCTs imply that the most probable time for a cell to divide is the current time.Once a cell has divided it is immediately able to divide again with high probability allowing cells proliferating under the exponentially distributed CCT assumption to reinforce their numbers.This is in direct contrast to cells with Erlang distributed CCTs (with the same mean but k 1) which, with high probability, will wait for a period of time before dividing.In short, the larger variance of the exponentially distributed CCT population allows it to grow more rapidly.number in each phase Fig. 3 The evolution of the of (a) the total numbers of cells and (b) the mean numbers of cells in each of k 4 stages.In panel (a) we plot the numerical (red line) and analytical (dashed black line) solutions for the total mean number of cells in the case k 4 and according to the naive (k 1) cell cycle model (analytical solution of equation ( 7) -blue line).In panel (b) analytically determined solutions (see equations ( 12)-( 15)) are plotted as dashed black lines and their numerical counterparts on top as solid coloured lines.The average CCT is C 10 arbitrary time units.The average period of each stage is equal (µi Cßk 2.5, λi 1ßµi 0.4 for i 1, . . ., k).
Under the assumption of identical transition rates, equation system (4) can be reduced to a closed equation for the mean number of cells in the k th stage and a set of k ¡ 1 ODEs which relate the number of cells in the other stages to Under the given initial conditions, a single cell in the first stage and no cells in any other stages, we can solve these equations to find where z is the first n th root of unity (Kendall, 1948).Although this expression looks complicated in some cases it is possible to express M j in a simple closed form.For example, when k 4 A comparison between these analytical solutions and their numerically solved counterparts demonstrates their veracity in Figure 3 (b).
By summing equation ( 11) over all values of j 1, . . ., k we can also find an expression for the total number of cells in a population: Although these formulae (equations ( 12)-( 16)) may be useful in specific cases where the closed form of the solution is readily accessible, their real utility is in shedding light on the asymptotic properties of the mean numbers of cells.
In the limit that t gets large for finite k the dominant term in the summation in equation ( 11) corresponds to the case r 0. Indeed for k 28 the real parts of the other elements in the summation are negative and hence these terms decay (Kendall, 1948).Thus we have where Summing equation ( 17) over all k stages leads to the asymptotic behaviour of the cell population as a whole: Equation ( 19) can be re-written as For all k 1, not only is the base of the exponent tßC less then e (since α k 1, for k 1), but the coefficient is less than unity (Kendall, 1948).This implies that, asymptotically, the expected total number of cells in a multi-stage model will always be less than the number expected in a single stage cell cycle model (which can be determined upon substituting k 1 in to ( 19)).
Note that in the limit as k , α k ln 2. Thus, as might have been expected for the deterministic model resulting from the limit k , the asymptotic population grows with base 2, doubling at regular intervals as the cells divide synchronously: Surprisingly though, the coefficient of 2 kßCt does not tend to unity in equation ( 21 ( Reversing the order of limits and taking the limit as k tends to infinity of equation ( 16) for finite t gives the limit for non-integer value of tßC, where ØxÙ gives the integer part of x (Kendall, 1948).For integer values of tßC the limit is corresponding to the algebraic mean of the limits of equation ( 23) as integers values are approached from the left and right hand sides.This "deterministic" doubling process is unsurprising since the waiting time distribution tends to a delta function in the k limit, implying that cell division is synchronous.

Cells are not distributed proportional to stage length
Returning to equations (4) under the assumption of identical rates of progression through the stages, we can derive corresponding equations for the mean proportion of cells in each stage, Mj M j ßM for j 1, . . ., k: At steady state we have the following recurrence relations for the mean proportion of cells in each stage In particular, this implies that M st j M st j¡1 for j 2 . . .k, so that, at steady state, as we progress through the stages we will have successively fewer cells in each stage on average (independent of the initial distribution of cells amongst different stages).By solving these recurrence formulae we can find exact expressions for the steady state proportions: In particular, note that M st That is to say there are twice as many cells in the first stage as the last stage at steady state when the number of stages is large.
These differences are potentially important for determining average CCTs experimentally.One popular method for determining cell cycle times is to label S-phase cells using two sequentially administered distinct DNA specific labels (Bokhari and Raza, 1992;Wimber and Quastler, 1963).
The administration of the labels is separated by a known time period.By counting cells labelled with one or both labels, and with reference to the known time period of separation, it is possible to calculate the mean duration of the S-phase.Once the proportion of cells in S-phase and the mean duration of S-phase have been determined it is also possible to calculate the mean cell cycle time for the population (Nowakowski et al., 1989).
The method outlined above implicitly makes the assumption that the number of cells in a particular phase of the cell cycle is proportional to the length of that phase.For the multi-stage model we have demonstrated that in the large time limit this is unequivocally not the case.Equation ( 11) can also be used to show that this phenomenon holds dynamically, although the mathematics is cumbersome.Instead we solve equations ( 4

The exponentially modified Erlang distribution
Although the identical-stage model, which gives rise to the Erlang distribution for CCTDs, is convenient from a mathematical perspective, it has been shown to have been outperformed by a number of other distributions (Golubev, 2010(Golubev, , 2016)).In particular, by considering 77 independent CCT data sets, Golubev (2016) has recently shown that one of the most appropriate distributions for representing CCTDs is the exponentially modified Gamma (EMG) distribution.For our pur- poses we will require that the shape parameter of the Gamma distribution is be integer valued so that the CCTD is actually an exponentially modified Erlang (EME).This will mean that we can continue to simulate CCTs using a series of exponentially distributed random variables (albeit one of them will have a different rate).Consequently this will allow us to continue to appropriately simulate processes in which cell division is important using the popular Gillespie algorithm.In order to generate EME distributed CCTs we modify our multi-stage cell cycle model as follows: Note that, in system (29), the rate of progression is identical through each of the initial k stages of cell cycle and that we have added an additional exponentially distributed stage at the end whose rate, λ 2 , is distinct from the rate, λ 1 , of the previous k stages.
We can ascertain the probability density function for the EME distribution, P EME ÔtÕ, by convolving the Erlang (P ER ÔtÕ) and exponential (P E ÔtÕ) distributions as follows, where λ 2 is the rate of the exponential distribution with which we are convolving and, as before, λ 1 is the rate of progression through each of the k identical exponentially distributed stages which comprise the Erlang distribution.We can simplify expression (30) to the following formulation where L λ 1 ¡λ 2 and Γ Ôk, LtÕ Lt z k¡1 e ¡z dz is the complementary incomplete gamma function.
We note that it is almost as simple to simulate this more general distribution in the Gillespie algorithm using a series of exponentially distributed stages as it is to simulate the distribution with constant rates of progression between stages.Indeed the simulation of any hypoexponential CCTD is straightforward in the Gillespie algorithm.However, the addition of extra parameters hampers their identifiability when fitting to experimental data and as such we only suggest using the Erlang or exponentially modified Erlang distributions in models of the CCTD.
In the following section we illustrate the importance of incorporating non-exponentially distributed CCTDs into stochastic simulations of cellular proliferation.For ease of understanding we concentrate purely on Erlang CCTDs, and note that the parameters are not based on fitted CCTDs but merely chosen for illustrative purposes.
Under the assumption of exponential CCTDs, Turner et al. (2009) write down and solve a simple probability master equation for the stem cell population.In particular, they consider the case in which r 1 r 3 which implies a positive net growth rate of the stem cell population.Under this assumption the mean number of cells in the stem cell population and variance can be shown to increase exponentially.Since the model is linear, by appealing to the central limit theorem, Turner et al. (2009) argue that, for large enough cell populations, the exact mean field equations given by dS dt will approximate the stochastic dynamics well.
In order to ensure a more realistic representation of the CCTD we introduce a multi-stage cell division process, as suggested above, so that the CCTDs are now Erlang distributed with the same mean ρ s , but with shape parameter k (and thus scale parameter µ 1 1ßÔρ s kÕ).As, before, at each division event, a choice about the type of division to occur is made with the same probabilities (r 1 , r 2 , r 3 ) as previously specified.Although we still get exponential increases in the mean and variance, the rate of increase is significantly decreased (see Figs. (c) Fig. 5 The evolution of the of the mean numbers, (a), variance, (b), and probability of the tumour having more than 1000 stem cells, (c), for cancer stem cells under model ( 32)-( 34) with varying numbers of exponentially distributed stages of the cell cycle, k 1, 2, 5, 10.The average CCT is ρs 1 and division probabilities are r1 0.2 and r3 0.15.Results are calculated from M 10, 000 repeat simulations.
Under this model with r 1 r 3 , if a tumour is not completely eradicated by treatment it is possible that it can return.It may, therefore be informative to know the probability that a tumour will reach a certain size by a particular time in order to plan appropriate follow-up therapeutic intervention.For example we may be interested to know the evolution of the probability that the tumour has reached 1000 stem cells in size (which we will denote p 1000 ÔtÕ) so that we might calculate the most appropriate time to initiate the follow-up intervention.In figure 5

Growth to confluence assays
Next we investigate the effect of incorporating a more realistic model of cell proliferation on the behaviour of a spatially extended individual-level model of cell migration and proliferation (Baker and Simpson, 2010).As such, we alter the mechanism of cellular proliferation from the original, exponentially distributed division times to our more realistic multi-stage Erlang distributed division times and observe the effect this has on the growth of the cell population.In order to achieve this we break the proliferation process into k stages, the passage through each of which has an exponentially distributed waiting time (as described above).As before, we chose the parameter of each stages' waiting time to ensure we have the same mean proliferation attempt time as in the original model.
In more detail, we consider a volume-exclusion process on a regular, square lattice in two dimensions with periodic boundary conditions.Each lattice site, of length h, can hold at most one In the multi-stage model, implementation (1) generally leads to less dense colonies than implementation (2) since cells do not attempt division as frequently.Under implementation (2) (see figure 6 (e) and (f)) a clear proliferating rim of (grey) cells can be seen with the bulk of cells being kept at stage k (black).Under implementation (1) every cell can be found in any stage of the cell cycle so it is hard to distinguish the proliferating rim (see figure 6 (b) and (c)).The difference between the two implementations however, is not due to aborted proliferation events in the bulk (away from the rim) but to the ability of cells at the proliferating rim to rapidly undergo a further division attempt after an aborted attempt under implementation (2).This suggests that the difference between the two implementations will only be apparent at high densities for which correlations have built up and significant numbers of division attempts are being aborted.
For low density systems, in which very few particles are adjacent, the mean cell division attempt times are almost the same for all values of k independent of the implementation ((1) or ( 2)).
However, the variance in the CCTDs for low density systems affects the rate of growth with larger values of k (less variance in the CCTD) generally leading to slower growth.This effect can be understood by considering equations ( 16) and ( 19) for a non-excluding population of cells, for which the finite time and asymptotic time behaviours, respectively, of cell populations with different values of k can be contrasted.
In figure 7  (f) Fig. 7 The influence of the number of proliferation stages, k, and proliferation rate Pp on the evolution of average cell density.Parameters, initial conditions and domain descriptions are as in figure 6. Panels (a)-(c) represent implementation (1) for Pp 0.05, 0.5, 1 respectively.Increasing the number of stages in the cell cycle causes growth of cell density to be retarded throughout the simulation.Panels (d)-(f) represent implementation (2) for Pp 0.05, 0.5, 1 respectively.Increasing k causes an initial retardation in growth followed by an acceleration as the effect of density correlations becomes prevalent.Note that figures in which Pp 1 are plotted on rescaled time axes for comparison (as described in the main text).
Under implementation (1) (see figure 7 (a)-(c)), even as cell-cell correlations build up, multistage cells still proliferate more slowly than single-stage cells, since the mean division attempt time remains the same for all values of k.The increased variance of cells with fewer stages results in faster population growth.
However, under the more realistic implementation (2) (see figure 7 (d)-(f)) cells with multistage cell cycles are able to re-attempt division after abortive division events more quickly than they otherwise could under the single-stage cell cycle model.Thus, effective average CCTs for cells with a multi-stage cell cycle at the proliferating rim of a cluster decreases in comparison to cells with a single-stage cell cycle.The faster the pairwise correlations build up, the more pronounced this effect becomes.With a very low proliferation rate (in comparison to fixed motility -see figure 7 (d)) cell movement is effective at breaking up correlations meaning that large clusters do not form and that we only see the effect of decreasing mean CCT for larger values of k at late (scaled) times, when the density is higher.Contrastingly, when proliferation is high in comparison to motility (see figure 7 (f)) clusters can form quickly preventing the bulk of cells from proliferating earlier and allowing the cells with multi-stage cell cycles to divide faster on average at the proliferating rim of these clusters than their single-stage counterparts.
It is also worth noting that the greater synchrony in cell division times for larger values of k (exemplified by equations ( 23)-( 24) for the limit of infinite k) is in evidence in the jagged nature of the yellow curves (corresponding to k 100) in all six subfigures.

Discussion
Currently, many stochastic models which incorporate cell proliferation employ the ubiquitous Gillespie stochastic simulation algorithm (Gillespie, 1976(Gillespie, , 1977)).Unfortunately, in its basic form the Gillespie algorithm represents all events as exponentially distributed.Cell cycle times are not exponentially distributed and can not, therefore be represented by a single reaction event in the Gillespie algorithm.Modelling cell cycle times as a single exponentially distributed event can lead to significant alterations in model behaviour in comparison to more appropriate CCTDs.Conse-of these correlation effects can be attributed to the environment in which the cells are proliferating.However, in NIH 3T3 cells a clear correlation has been observed between daughter cells of a given mitotic event compared to more distant relatives; implying a heritable predisposition (Mort et al., 2014).Therefore, one obvious extension to this work would be to incorporate the effects of correlations in cell and phase times to better reflect the biological heterogeneity of a given system.

Appendix A -Materials and Methods
In order to determine cell cycle times in cell culture, NIH 3T3 Flp-In cells (Mort et al., 2014) were seeded at various densities in phenol-red free dulbeccos modified eagle medium (DMEM) containing 10% fetal calf serum, 1% Penicillin/Streptomycin and 100 µgßml Hygromycin B on glass bottomed 24-well culture plates (Greiner bio-one, UK).The next day, time-lapse imaging was performed on subconfluent wells with a 20¢ objective using a Nikon A1R inverted confocal microscope in a heated chamber supplied with 5% CO 2 in air.The time elapsed between mitotic events was measured using the Fiji (Schindelin et al., 2012) distribution of ImageJ -an open source image analysis package based on NIH Image (Schneider et al., 2012).

Appendix B -Psuedocode for the multi-stage model of the cell cycle
One of the original (and most popular) implementations of the mathematically exact SSA is known as the direct method (Gillespie, 1977).Here we present pseudocode for the direct method implementation of the simple multi-stage model of the cell cycle (corresponding to Erlang distributed CCTs) in a well mixed context.
Let XÔtÕ be a vector of length k which contains the number of cells in each stage at time t.A time interval τ , until the next stage advancement event, is generated.Along with it, an index j, is chosen which determines from which stage a cell will advance from time t τ .The changes in the numbers of cells caused by the stage advancement are implemented, the propensity functions (progression rate, λ, multiplied by the number of cells in each stage) are altered accordingly and (see Fig 2 for one such example).

Fig. 2
Fig.2The agreement between experimentally determined cell cycle times in NIH 3T3 mouse embryonic fibroblasts grown in vitro (grey histograms) and (a) the Erlang distribution (green curve) (b) the exponentially modified Erlang distribution (blue curve).

Figure 2
Figure 2 demonstrates the improved agreement between the Erlang and exponentially modified ) as might have been expected.Thus the total population grows like lim k lim t M 0.721 ¤ 2 ktßC .
) numerically.Numerical solution also allows us to investigate the more general hypoexponential CCTD model (2) for which no analytical solutions are available.Figures 4 (a) and (b) display the evolution of the mean numbers and proportions (respectively) of cells in each stage for equal stage progression rates, λ 1 , and figures 4 (c) and (d) display the equivalent for unequal progression rates.The number/proportion of cells in each stage is not proportional to the mean duration of the stage, µ i , either at steady state (compare actual steady state proportions with the stars representing the normalised mean stage durations, µ i ß µ i , in 4 (b) and 4 (d)) or dynamically.

Fig. 4
Fig. 4 The evolution of the of the mean numbers ((a),(c)) and proportions ((b),(d)) of cells in each of k 5 stages.In all cases the average CCT is C 10 arbitrary time units.Panels (a) and (b) represent the scenario in which rates are chosen so that the average period of each stage is equal (µi Cßk 2, λi 1ßµi 0.5 for i 1, . . ., k), (c) and (d) the scenario in which the mean stage durations, µi are chosen by partitioning the total CCT uniformly at random into k 5 parts.Rates are given by λi 1ßµi for i 1, . . ., k. Stars at t 30 on (b) and (d) indicate the expected proportion of cells at steady state if numbers of cells in each stage were proportional to cell stage duration.For the equal progression rates model all stars overlap in (b).
5 (a) and (b)).Crucially this means that the deterministic mean-field model derived from the original process will significantly overestimate the number of cancer stem cells in the tumour which could have significant therapeutic implications.
(c)  we plot the evolution of p 1000 ÔtÕ over time.It is clear, by t 100, that the probability of the tumour having grown to 1000 stem cells, p 1000 varies significantly depending on the value of k used in the model despite the cells having the same mean CCT 5 .The effects of varying CCTD can clearly be seen to influence the model outcome even in this relatively straightforward linear model of cellular proliferation.In more complex models, in which other species depend in a non-linear manner on the number of cells, the effects will no doubt be further exacerbated.The potential for therapeutic interventions to be based on stochastic mathematical models of cellular proliferation further emphasises the importance of modelling the CCTD correctly.

Fig. 6
Fig. 6 The influence of the number of proliferation stages, k, and the proliferation abortion mechanism on the spatial coverage of cells populating the domain at t 10.Cells in different stages are represented by different shades of grey.Darker shading corresponds to later stages and white indicates empty sites.Parameters are Pm 1, Pp 1 with Lx Ly 100.Initial seeding density is 1% (i.e. 100 cells).Panels (a)-(c) represent implementation (1) for k 1, 10, 100 respectively.Increasing the number of stages in the cell cycle causes a decrease in terminal cell density in this scenario.Panels (d)-(f) represent implementation (2) for k 1, 10, 100 respectively.Increasing k causes an increase in the terminal cell density in this scenario.
we contrast the evolution of the spatially-averaged density for three values of k 1, 10, 100 and three proliferation rates, P p 0.05, 0.5, 1, under implementation (1) ((a)-(c)