Royal Society Publishing

An automaton model for the cell cycle

Atilla Altinok, Didier Gonze, Francis Lévi, Albert Goldbeter

Abstract

We consider an automaton model that progresses spontaneously through the four successive phases of the cell cycle: G1, S (DNA replication), G2 and M (mitosis). Each phase is characterized by a mean duration τ and a variability V. As soon as the prescribed duration of a given phase has passed, the transition to the next phase of the cell cycle occurs. The time at which the transition takes place varies in a random manner according to a distribution of durations of the cell cycle phases. Upon completion of the M phase, the cell divides into two cells, which immediately enter a new cycle in G1. The duration of each phase is reinitialized for the two newborn cells. At each time step in any phase of the cycle, the cell has a certain probability to be marked for exiting the cycle and dying at the nearest G1/S or G2/M transition. To allow for homeostasis, which corresponds to maintenance of the total cell number, we assume that cell death counterbalances cell replication at mitosis. In studying the dynamics of this automaton model, we examine the effect of factors such as the mean durations of the cell cycle phases and their variability, the type of distribution of the durations, the number of cells, the regulation of the cell population size and the independence of steady-state proportions of cells in each phase with respect to initial conditions. We apply the stochastic automaton model for the cell cycle to the progressive desynchronization of cell populations and to their entrainment by the circadian clock. A simple deterministic model leads to the same steady-state proportions of cells in the four phases of the cell cycle.

1. Introduction

Remarkable advances were made during the last 20 years in elucidating the nature of the regulatory network that drives the cell division cycle. Building on experimental studies of the early cell cycles in amphibian embryos and of the yeast cell cycle, investigations of the mammalian cell cycle have shown that it is driven by a complex network of cyclin-dependent kinases (Cdks) [1]. Kinetic models for the cell cycle were first proposed for amphibian embryos [24] and yeast [5], and later for the mammalian cell cycle [68]. The early models for amphibian embryos typically counted a few variables, while models for the yeast and mammalian cell cycle incorporate a much larger number of variables. Thus, a model recently proposed for the Cdk network driving the mammalian cell cycle counts as many as 39 variables [6], which increases to 44 upon incorporation of a DNA replication checkpoint.

Several ways exist for reducing the complexity of models for the cell cycle. One consists in the construction of simpler models in which the number of variables is reduced by relinquishing a variety of biochemical details, which is necessary to achieve such a simplification. Such an approach was followed, for example, by Csikász-Nagy et al. [9] and by Gérard & Goldbeter [10], who studied the dynamics of a five-variable model that can be seen as a skeleton version of their much more detailed model for the Cdk network [6]. In an alternative approach, not focusing on the biochemical intricacies of intracellular dynamics, cell populations of various ages are described by continuous, partial differential equations [11].

Yet another approach rests on the use of a discrete automaton model for the cell cycle. The goal of this paper is to present such an automaton model for the cell cycle and to explore its dynamical properties. Our aim is to provide a simple description, free of biochemical details, in terms of an automaton switching between the successive phases of the cell cycle. Keeping in mind that the transitions between the various phases of the cell cycle have a random nature [1214], we incorporate into the cell cycle automaton model the stochastic aspect of the transitions between successive cell cycle phases. We show that this stochastic automaton model for the cell cycle can be used to study the progressive desynchronization of cell populations and to account for their entrainment by an external signal such as the circadian clock. Another application pertains to cancer chronotherapy [15,16]. Probing the effect of anti-cancer drugs by means of modelling and numerical simulations requires a model for the cell cycle. To do this, one may use a detailed biochemical model of the cell cycle [6], or a version containing a reduced number of variables [10]. Such models are useful, for example, in addressing the effect of mutations, or of inhibitors of enzymes such as Cdks, which control the transitions between cell cycle phases. However, when the distribution of cells between the various cell cycle phases becomes the predominant issue, biochemical details do not play such a major role, and the automaton model provides a useful, versatile tool for describing the dynamics of a population of proliferating cells. Thus, we previously used the cell cycle automaton model to study the response of a population of proliferating cells to the circadian administration of an anti-cancer drug that acts on a particular phase of the cell cycle [1719].

Here, we analyse in detail the dynamics of the automaton model and its general properties. We examine how its dynamics is affected by factors such as the mean durations of the cell cycle phases and their variability, the type of distribution of the durations, the number of cells, the regulation of the cell population size and the initial conditions. We show how the distribution of cells can be related to the DNA content. Finally, we present a simple deterministic model that leads to the same steady-state proportions of cells in the four phases of the cell cycle.

2. An automaton model for the cell cycle: assumptions and algorithm

The automaton model for the cell cycle is shown in figure 1a. It is based on the following assumptions:

  • — The cell cycle consists of four successive phases along which the cell progresses: G1, S (DNA replication), G2 and M (mitosis). The possibility of temporary arrest in a G0 phase is not considered here but the automaton can easily be extended to include such an additional phase.

  • — Each phase is characterized by a mean duration τ and a variability V. Typical mean durations of the cell cycle phases are shown in figure 1b. As soon as the prescribed duration of a given phase has passed, the transition to the next phase of the cell cycle occurs. The time at which the transition takes place varies in a random manner according to a distribution of durations of the cell cycle phases. In the case of a uniform probability distribution, the duration varies in the interval [τ (1 − V), τ (1 + V)]. We also consider the case of a lognormal distribution and compare it with the uniform distribution (see appendix A.2).

  • — Upon completion of the M phase, the cell divides into two cells, which immediately enter a new cycle in G1. The duration of each phase is reinitialized for the two newborn cells.

  • — At each time step in any phase of the cycle, the cell has a certain probability of being marked for exiting the cycle and dying at the nearest G1/S or G2/M transition. To allow for homeostasis, which corresponds to the maintenance of the total cell number, we further assume that cell death counterbalances cell replication at mitosis. The cell death probability must be of the order of 50 per cent over one cycle. When the probability of cell death is slightly smaller or larger than the value yielding homeostasis, the total number of cells increases or decreases in time, respectively. To stabilize the population around an arbitrary value, an option to allow regulation of the death probability by the population size is also included (see §3.2 below).

Figure 1.

(a) Scheme of the automaton model for the cell cycle. Each cell switches sequentially between the phases G1, S, G2 and M, after which the cell divides and the two daughter cells enter a new G1 phase. The duration of each phase is determined stochastically (according to a uniform or lognormal distribution). When the time of the current phase is elapsed, the cell switches to the next phase. At any time, each cell can be marked for death with a given probability, but the cells marked for death exit the cell cycle only at the next G1/S or G2/M transitions. (b) Default mean duration of each phase of the cell cycle considered in numerical simulations of the automaton model.

The algorithm corresponding to these assumptions is shown in figure 2. The code of the computer program based on this algorithm and used for numerical simulations of the cell cycle automaton model is given in the electronic supplementary material.

Figure 2.

Algorithm of the stochastic automaton model for the cell cycle shown in figure 1 (see also text for a detailed description). The corresponding computer code (in Matlab) is provided in the electronic supplementary material.

Although the model does not incorporate any explicit biochemical detail, we note, from a biological viewpoint, that the transitions between successive cell cycle phases are controlled by changes in the activity of various Cdks. Thus, the cyclin D/Cdk4–6 and cyclin E/Cdk2 complexes promote progression in G1 and elicit the G1/S transition, respectively; the cyclin A/Cdk2 complex ensures progression in S and the transition S/G2, while the activity of the cyclin B/Cdk1 complex brings about the G2/M transition. A link, therefore, exists between the automaton model and the kinetic models describing the changes in the activity of Cdks. In the automaton model, a cell progresses in a given cell phase until it switches to the next phase. In models that incorporate biochemical details, the transition from one phase to the next is associated with the transient activation of the corresponding Cdk. The automaton model encompasses a stochastic aspect that has yet to be incorporated into the detailed models for the Cdk network. Stochastic simulations of the latter models, which often contain a large number of variables, are more cumbersome than those carried out with the automaton model. The automaton model, therefore, provides a tool of choice for the study of desynchronization in a cell population and its response to phase-specific drugs that interfere with the dynamics of the cell cycle.

The automaton model for the cell cycle combines a Markovian process—cell death can occur at any time, independently of the history of the cell—with a non-Markovian process—the duration of each phase is pre-defined, according to a probability distribution. In spite of its mixed nature, such a model captures basic properties of the cell cycle, and numerical simulations produce realistic and interpretable behaviours. In an alternative approach, to be considered in future work, instead of specifying pre-defined phase durations, the transitions between the successive phases of the cell cycle are themselves stochastic, making the model consistently Markovian. We show below that the automaton model based on the above assumptions captures key features of the cell cycle and leads to specific predictions as to the effect of the variability in the duration of cell cycle phases.

3. Dynamics of the cell cycle automaton: from single cells to cell populations

3.1. Dynamics of a single cell

The dynamics of the cellular automaton for a single cell is illustrated in figure 3a. This figure shows the time spent by the cell in each phase of the cell cycle during successive cycles. The duration of each phase is chosen probabilistically from uniform distributions characterized by phase-specific means (τG1 = 9 h, τS = 11 h, τG2 = 1 h and τM = 1 h; figure 1b) and a variability V. We consider in the following that this variability is the same for each phase, but this does not need at all to be the case. The basal value of the mean duration of the cell cycle is of the order of those typically found for mammalian cells, which are often close to 24 h, although the duration of the cell cycle can markedly vary for different cell types. Thus, a study of squamous cell carcinoma of rat trachea reported the values τG1 = 15 h, τS = 9 h, τG2 = 3 h and τM = 1 h for a cell cycle of 28 h [20], while a recent study using phase-specific fluorescence biosensors in various human cell lines such as NIH3T3 and Hela cells gave values of the order of τG1 = 400–500 min, τS = 500 min, τG2 = 100 min and τM = 60 min for a cell cycle duration close to 19 h [21]. For a given cell type, the mean durations of some cell cycle phases may vary considerably [21]. In a given organism, these durations may also vary as a function of nutritional conditions, as shown in rat (e.g. [22]), or in the course of cell differentiation [23]. Together with their variability V, the mean durations of the various phases represent key parameters of the automaton model that can readily be modified. We have previously shown that the response to circadian administration of anti-cancer drugs markedly depends on the total duration of the cell cycle and the mean durations of its phases [19].

Figure 3.

(a) Dynamics of a single automaton undergoing repetitive cell cycles. (b) Distribution of the total cell cycle duration over 10 000 cycles. The duration of each cell phase is taken from uniform distributions with the means τG1 = 9 h, τS = 11 h, τG2 = 1 h and τM = 1 h (the total duration of the cell cycle is equal to 22 h), and a variability V = 20%.

The distribution of the total cell duration is centred around τc =τG1 +τS + τG2+ τM. In the case considered in figure 1b, for V = 20%, the cell cycle duration ranges from 18 to 26 h (figure 3b). Its distribution is not uniform: in the case considered in figure 1b, where the M and G2 phases are brief compared with G1 and S, the distribution of τc is close to a triangle, obtained from the sum of the uniform distributions taken for the durations of the G1 and S phases.

3.2. Ensuring the homeostasis of the cell population

At the end of one cell cycle, two cells emerge from one cell. In the absence of cell death, the number of cells would rapidly explode. To ensure homeostasis, 50 per cent of cells must exit the cell cycle during one cycle. Thus, each cell should have a probability of 0.5 of dying during each cycle. When this probability is higher, the cell population decreases exponentially. Conversely, a smaller value for the probability of exiting the cell cycle corresponds to an exponential increase of the population.

In the automaton model, at each time step, the cells have a certain probability of exiting the proliferating pool at the next G1/S or G2/M transition (they can be marked for death at each time step, but they are removed from the cell cycle only at the G1/S or G2/M transition). This assumption reflects the role of the molecular checkpoints in eukaryotic cells.

Homeostasis can be obtained in at least two ways. First, by adjusting the death probability P0, i.e. the probability of the cell being marked for death at each time step. Upon appropriate calibration of this parameter (P0 = Ph), stabilization of the cell number in the population can be obtained (figure 4a). The exact calculation of Ph is given in appendix A.1. If this probability is even slightly higher (P0 > Ph) or lower (P0 < Ph), then the cell population will, respectively, decrease (figure 4b) or increase (figure 4c) exponentially. When P0 = Ph, the population stabilizes around the initial number of cells N0.

Figure 4.

Homeostasis and the probability of cell death. (ac) The probability of death is not regulated by the size of the population and the number of cells either stabilizes at its (a) initial value, (b) decreases or (c) increases, depending on the relative value of P0 with respect to Ph (see appendix A.1). (df) The probability of cell death is modulated by the cell population size according to equation (3.11). In (a, d), P0 = 0.525 × 10−3 = Ph. In (b,e) P0 = 0.55 × 10−3 > Ph. In (c,f) P0 = 0.50 × 10−3 < Ph. In (df) ks = 0.001 min−1, NS = 104. At t = 0, the total number of cells is 104 and all cells start in the G1 phase. The total cell cycle duration is equal to 22 h (τG1 = 9 h, τS = 11 h, τG2 = 1 h and τM = 1 h), and the variability V = 20%.

Alternatively, to stabilize the cell population around a given value, e.g. 104 cells, which is the typical number of cells considered in numerical simulations, a regulation of the cell population may be introduced via a logistic equation that incorporates the total number of cells, N, as well as the target number of cells, Ns, around which the population should stabilize,Embedded Image 3.1where P0 = rΔt is the probability of dying at each time step Δt, r being the death rate (in min−1). In all the simulations presented in the paper, Δt = 1 min and thus the probability P0 is defined accordingly; moreover, ks will be taken to be between 10−3 and 10−4 min−1 and Ns = 104 cells. As shown below, the value of parameter ks affects the rate of stabilization of the population around Ns, but it must be sufficiently reduced so that the logistic term remains small with respect to P0. When the population is above Ns, the propensity to quit the cell cycle increases to stabilize the population around the steady-state value Ns. Conversely, this propensity decreases when the population is under Ns.

As illustrated in figure 4df, using the logistic equation (3.1) ensures homeostasis without requiring any fine-tuning of the death probability P0. Furthermore, the population converges to Ns regardless of the initial value N0.

3.3. Effect of the number of cells

After dealing with the dynamics of a single cell, we may now apply the cell cycle automaton model to the dynamics of a cell population, under the conditions of homeostasis outlined in §3.2. We shall focus on the case where the cells proceed in the cell cycle independently from one another. We first determine the effect of the initial number of cells, when homeostasis is achieved by setting P0 = Ph (the number of cells then remains equal to their initial number; see §3.2). The dynamics of the automaton for cell populations of various sizes is shown in figure 5 for cell numbers ranging from 10 (figure 5a) to 10 000 (figure 5d). All cells start at the beginning of phase G1. The variability V of the duration of each phase is 20 per cent. Owing to the random variations in the duration of the phases, cells do not switch to the next phase all at the same time. This results in irregular oscillations in the fractions of cells in each phase. We shall see in the next section that these oscillations possess a transient nature: the oscillations dampen and the system eventually reaches a steady-state distribution of cells in each phase. As the population size increases, the oscillations become more regular, since random fluctuations progressively average out.

Figure 5.

Dynamics of the cell cycle automaton for increasing values of the initial number of cells in the population: (a) N0 = 10, (b) N0 = 100, (c) N0 = 1000 and (d) N0 = 10 000. In all cases, the total cell cycle duration is equal to 22 h and the mean durations of the phases are as in figure 4. Moreover, V = 10%, P0 = 0.525 × 10−3 (homeostasis is achieved without logistic regulation of the probability of cell death).

3.4. Effect of variability of durations of cell cycle phases

The variability in the duration of the cell cycle phases is responsible for progressive cell desynchronization. In the absence of any variability (V = 0), if we assume that the duration of each phase is the same for all cells and that all cells start at the same point of the cell cycle, e.g. at the beginning of G1, a sequence of square waves that bring the cells synchronously through G1, S, G2 and M, and back into G1, occurs (figure 6a). The square waves will continue unabated over time. However, as soon as some degree of variability of the cell cycle phase durations is introduced, these square waves transform into oscillations through the cell cycle phases, the amplitude of which diminishes as the variability increases, as shown in figure 6bd for V = 5%, 10% and 20%, respectively. In the long term, these oscillations dampen as the system settles into a steady-state distribution of cell cycle phases: the cells are fully desynchronized and have forgotten that they initially started from the same point of the cell cycle. The damping of the oscillations becomes faster as the variability increases, as shown over a longer time course in figure 6e,f for V = 20% and 50%.

Figure 6.

Effect of the variability of the durations of the cell cycle phases: (a) V = 0%, (b) V = 5%, (c) V = 10%, (d) V = 20%, (e) V = 20% (long time run) and (f) V = 50% (long time run). The mean durations of the cell cycle phases are as in figure 4; P0 = 0.525 × 10−3 (homeostasis is achieved without logistic regulation of the probability of cell death). The initial number of cells is N0 = 10 000 cells.

3.5. Steady-state distribution of cell cycle phases: independence of initial conditions

The initial distribution of the cells in the different phases may affect the transient oscillations but not the steady-state distribution of cells in the different phases. In figure 7, we show the time evolution of the fraction of cells in each phase, for the model without (figure 7a,b) and with (figure 7cf) logistic regulation of the cell number, starting from different initial conditions: all cells in G1 (figure 7a,c), or 25 per cent of the cells in each phase (figure 7b,d). We see that the same steady-state distribution of cells in the four phases of the cell cycle is established regardless of initial conditions. The same result is obtained when simulations of the cell cycle automaton model start with a single cell in G1, but the evolution is slower when the value of ks in equation (3.1) is reduced (compare figure 7e with figure 7f).

Figure 7.

Effect of initial conditions: independence of the steady-state proportions of cells in the four phases of the cell cycle. Total cycle length and durations of the cell cycle phases are as in figure 4, with V = 50%, P0 = 0.525 × 10−3. In (a,b) homeostasis is achieved without logistic regulation of cell death probability: P0 = 0.525 × 10−3. In (cf) cell death probability is regulated according to the logistic equation (3.1), with P0 = 0.525 × 10−3, ks = 0.001 min−1 (ce) or 0.0002 min−1 (f), NS = 104. (a,c) All cells are initially in G1. (b,d) Initially, 25% of cells are in each phase. (e,f) Initially, there is a single cell in G1. In all cases the automaton evolves towards the same steady-state proportions of cells in each phase.

4. Predicting steady-state proportions of cells in the various phases of the cell cycle

It is useful to develop an alternative approach to predict the steady-state distributions of cells in each phase. To this end we use the following deterministic model, which contains four variables that represent the amounts of cells in each phase:Embedded Image 4.1

Parameters τM, τG1, τS and τG2 denote, as before, the duration of the corresponding cell cycle phases. The probabilities of a cell not dying in a given phase are given byEmbedded Image 4.2where Ps is the survival rate. Consistent with the conditions for homeostasis in the automaton model (see §3.2), we haveEmbedded Image 4.3or, when using the logistic regulation,Embedded Image 4.4

Equations (4.1) yield a single steady-state distribution of cells between the four phases. This steady-state distribution is in agreement with the predictions of the automaton model (figure 8). Thus, both in the deterministic model and in the stochastic automaton model, the proportion of cells in each phase only depends on the respective durations of the cell cycle phases.

Figure 8.

Effect of the mean cell cycle phase durations on the steady-state proportions of cells in each phase. The total cell cycle duration is equal to 22 h in (a), with τG1 = 9 h, τS = 11 h, τG2 = 1 h and τM = 1 h. In (b) the total cell cycle duration is still equal to 22 h, but τG1 = 14 h, τS = 6 h, τG2 = 1 h and τM = 1 h. Homeostasis is achieved without regulation of the probability of cell death (P0 = Ph = 0.525 × 10−3), N0 = 1000 and V = 20%. The steady-state proportions of cells in each phase are the same as those predicted by the deterministic model in §4. When setting dG1/dt = dS/dt = dG2/dt = dM/dt = 0 in equations (4.1), we obtain the following proportions at steady state: G1 = 0.49, S = 0.45, G2 = 0.029 and M = 0.028 for the mean phase durations of (a), and G1 = 0.74, S = 0.20, G2 = 0.028 and M = 0.027 for (b).

5. Entrainment of the cell cycle automaton by the circadian clock

Experiments indicate that the cell cycle in many tissues can be entrained by the circadian clock [2426]. Molecular links between these two oscillators have indeed been identified [16]. Thus, the kinase Wee1 that controls the G2/M transition is induced by the circadian clock [27]. It is useful to investigate how the automaton model for the cell cycle can account for its entrainment by the circadian clock. Entrainment can be included in the automaton model by considering that the protein Wee1 undergoes circadian variation owing to induction of the wee1 gene by the circadian clock (figure 9). Wee1 is a kinase that phosphorylates and thereby inactivates the Cdk, Cdk1 that controls the G2/M transition and, consequently, the onset of mitosis.

Figure 9.

Scheme of the model showing the control of the cell cycle by the circadian clock via Wee1 and Cdk1, which, respectively, inhibits and induces the G2/M transition.

Because the cell cycle automaton model can be applied to probe the physiological effect of the timing of drug administration, let us focus on the case of the human circadian clock when investigating its coupling to the cell cycle. In general, humans adhere to a routine of 16 h of diurnal activity alternating with 8 h of nocturnal sleep; thus, in the case of human modelling, we will consider a 16 L : 8 D cycle (16 h of light, from 8.00 to 24.00 h, followed by 8 h of darkness, from 24.00 to 8.00 h). The rise in Wee1 should occur at the end of the activity phase. The decline in Wee1 activity is followed by a rise in the activity of the kinase Cdk1, which enhances the probability of transition to the M phase. We shall consider that the rise in Wee1 is immediately followed by a rise in Cdk1 (figures 9 and 10).

Figure 10.

Entrainment of the cell cycle automaton by the circadian clock via the circadian variation of Wee1 and Cdk1, which varies in a semi-sinusoidal manner. Wee1 varies in a semi-sinusoidal manner with a peak at 22.00 h so as to increase between 14.00 and 22.00 h and decrease from 22.00 to 4.00 h; it remains equal to zero from 4.00 to 14.00 h. Cdk1 varies in a similar manner, with a lag of 6 h, so as to increase between 22.00 and 4.00 h and decrease from 4.00 h to 22.00 h; it remains equal to zero from 10.00 to 22.00 h. The duration of the cell cycle prior to entrainment is equal to 22 h (τG1 = 9 h, τS = 11 h, τG2 = 1 h and τM = 1h) in (a), and 26 h (τG1 = 10 h, τS = 12 h, τG2 = 2 h and τM = 2 h) in (b). For both cases, V = 20% and homeostasis is achieved by logistic regulation of the probability of cell death, with P0= Ph = 0.000525, Ns = 10 000 and ks = 0.001.

In the cell cycle, we consider that the probability of transition from G2 to M, at the end of the G2 phase, decreases as Wee1 rises, according to equation (5.1),Embedded Image 5.1Conversely, we assume that the probability of premature transition from G2 to M (i.e. before the end of G2, the duration of which was set when the automaton entered in G2) increases with the activity of Cdk1 according to equation (5.2),Embedded Image 5.2The probability is first determined with respect to Cdk1; if the G2/M transition has not occurred, the cell progresses in G2. Only at the end of G2 is the probability of transition to M determined as a function of Wee1. We assume here that Wee1 varies in a semi-sinusoidal manner with a peak at 22.00 h and that Cdk1 varies in a similar manner, with a lag of 6 h, so as to reach a peak at 4.00 h. Other patterns mimicking the situation occurring in other species such as the mouse may easily be implemented.

Upon entrainment by the circadian clock, cells become more synchronized than in the absence of entrainment. In the case considered in figure 10a, the period changes from 22 to 24 h, which corresponds to the period of the external LD cycle. In contrast to the progressive dampening of the oscillations in the absence of entrainment (figure 6e,f), oscillations become sustained when the cell cycle automaton is driven by the circadian clock. The model shows that the cell cycle can also be entrained to a circadian period when the cell cycle duration exceeds 24 h. Thus, a 26 h-long cell cycle can switch to a 24 h duration when coupled to the circadian rhythm in Wee1 and Cdk1 (figure 10b). Here again the entrainment ensures that the oscillations in the fractions of cells in the four phases of the cell cycle become sustained.

6. Discussion

We studied, in this work, the general properties of an automaton model for the cell cycle. The model is based on the stochastic transitions between successive phases M, G1, S and G2 of the cell cycle. The duration of each phase varies randomly from cell to cell around a mean value. When applied to a population of initially synchronized cells, i.e. a population of automata with different durations of cell cycle phases, such a variability leads to transient oscillations in the evolution of the fractions of cells in each phase until, finally, a steady-state distribution of cells between the four phases is established (figure 6). By specifying that the probability of cell death is modulated by cell density, the automaton model can be tuned to ensure that the cell population stabilizes at a prescribed level, thus allowing for homeostasis.

We illustrated the effect of the number of cellular automata on the dynamics of a population. While fluctuations in the fractions of cells in the four phases are conspicuous at low cell numbers, they are significantly attenuated when the population count is only a few hundreds of cells (figure 5). The amplitude of the transient oscillations and the time required for reaching the steady-state distribution both decrease when the variability in the durations of the cell cycle phases increases (figure 6). We further showed that the final proportions of cells in each phase at steady state are independent of the initial conditions: the same distribution is indeed reached when cells all start in a given phase, or when the population starts with similar amounts of cells in each phase, or even with a single cell (figure 7). A deterministic model predicts the same steady-state proportions of cells in each phase as those obtained by means of the automaton; these proportions vary with the mean durations of the cell cycle phases (figure 8). Numerical simulations were performed throughout this paper for the case of uniform distributions of the durations of cell cycle phases. Similar results are obtained when using other descriptions for the random variations, such as the lognormal distribution (see appendix A.2 and the electronic supplementary material, figure S1).

The distribution of cells between the various phases of the cell cycle can be characterized experimentally by means of the fluorescence-activated cell sorting (FACS) methodology, according to their DNA content in the different phases. The automaton model can be used to generate ‘virtual’ FACS profiles showing the dynamical changes of the distribution of cells between the various cell cycle phases, assuming that the quantity of DNA increases linearly from the beginning of S with the time spent in this phase, until it doubles at the end of S (see electronic supplementary material, figure S2).

In many organisms, including humans, the cell cycle is coupled to the circadian clock. In the presence of such a coupling, the automaton model for the cell cycle shows that oscillations in the fractions of cells in the four phases become sustained [17,18]. This result accounts for the observation that cells of many tissues in humans display a circadian peak in the time course of the fractions of cells in the four phases [24,26]. The cell cycle automaton can be entrained by the circadian clock, whether its autonomous duration is less or greater than 24 h (figure 10). Other theoretical studies have addressed the dynamical consequences of interactions between the cell cycle and the circadian clock [28].

The term ‘cellular automaton’ is often used to describe the spatio-temporal evolution of chemical or biological systems that are capable of switching sequentially between several discrete states [29,30]. In regard to cell proliferation and morphogenesis, cellular automata have been used to model tumour growth [31,32] and pattern formation [33]. One typical application of cellular automata pertains to excitable systems (e.g. neurons or muscle cells), which can evolve from a rest state to an excited state, then to a refractory state, before returning to the rest state [34]. Spatially coupled automata can account for the propagation of waves in excitable media. Here, as in a follicular automaton model for hair cycles [35,36], we consider spatially independent automata. Thus we assume that, within a population, the dynamics of a cycling cell will not be influenced by the state of neighbouring cells. In contrast, in excitable systems, a cell in the rest state can be triggered to switch to the excited state when a neighbouring cell is excited. Keeping these differences in mind, we note that the model for the cell cycle can be viewed as a true ‘cellular automaton’, if only because it pertains to the cell cycle!

When considering homeostasis of the cell population in §3.2, we introduced implicitly a form of cell-to-cell communication since the dynamics of the automaton is controlled by the size of the population through the parameter measuring cell death. The automaton model could readily be extended to the case of a spatial array of cells in which the dynamics of a given cell would be affected by its neighbours, e.g. through contact inhibition. By taking into account explicitly such interactions between cells, the density effect would emerge naturally from the cell cycle dynamics. Such an approach could prove useful in studying the dynamics of cell proliferation in tissues or tumours.

The present model can serve various applications, including the study of the effect of anti-cancer drugs that kill cells at specific phases of the cell cycle. It was previously used to account for the differential cytotoxicity of 5-fluorouracil (5-FU), an anti-cancer drug that kills cells in the S phase, when it is administered according to a circadian schedule that peaks at 4.00 h (minimum toxicity) or 16.00 h (maximum toxicity) [37,38]. The automaton model shows that, upon entrainment by the circadian clock through Wee1 and Cdk1, the peak of the fraction of cells in S phase occurs around 16.00 h, so that the maximum number of cells will be killed when the peak of 5-FU delivery occurs at that time [17,18]. The opposite result is obtained when 5-FU peaks around 4.00 h, since the fraction of cells in S phase then goes through a trough (figure 10).

If 5-FU cytotoxicity is minimal around 4.00 h, how can one explain that the circadian administration of the drug peaking at that time is nevertheless more effective against tumour cells? We have used the cell cycle automaton model to study the effect of 5-FU on two cell populations, one of healthy cells and one of tumour cells, that differ by the cell cycle length, the mean durations of the cell cycle phases, the variability of these durations and the entrainment of the cell cycle by the circadian clock. Each of these differences can enhance the differential cytotoxicity of 5-FU towards the two cell populations when a single circadian pattern of 5-FU peaking at 4.00 h is considered. All these effects are additive, so that the differential cytotoxicity can be maximized when they are combined [18,19]. The automaton model thus helps us to uncover plausible mechanisms for the observation that a given circadian pattern of anti-cancer drug delivery can be, at the same time, least cytotoxic to healthy cells and highly toxic for tumour cells.

Acknowledgements

This work was supported by grant no. 3.4607.99 from the Fonds de la Recherche Scientifique Médicale (F.R.S.M., Belgium), by the Belgian Federal Science Policy Office (IAP P6/25 ‘BioMaGNet’: ‘Bioinformatics and Modelling: From Genomes to Networks’) and by the European Union through the Network of Excellence BioSim, contract no. LSHB-CT-2004-005137.

Footnotes

  • Received October 10, 2010.
  • Accepted November 1, 2010.

Appendix

A.1. Homeostasis of the cell population: calculation of Ph

To achieve homeostasis, the probability of a cell dying during one cycle should be equal to 0.5. The probability of dying during one cycle should thus be equal to the probability of surviving during the whole cycle,Embedded Image A 1where Δt is the sampling time (=1 min in all the simulations), T is the duration (in min) of a complete cell cycle and P0 = r Δt is the probability of dying at each time step, r being the death rate.

From this equation, we can calculate P0,Embedded Image A 2

Thus, for a cell cycle of total duration T = 22 h (=1320 min), the time step being fixed to Δt = 1 min, we findEmbedded Image A 3Note that this value is close to the value P0 = 0.0005380 previously found in an empirical manner through numerical simulations by Altinok et al. [17,18].

A.2. Uniform versus lognormal distribution of phase durations

All the results presented in this paper have been obtained for the case where the durations of the various cell cycle phases obey a probability distribution centred around the mean duration τ, with a range of variation extending uniformly from τ (1 − V) to τ (1 + V). In the program, we also implemented an option to generate phase duration according to a lognormal (rather than a uniform) distribution. The lognormal distribution obeys the following equation:Embedded Image A 4where μ denotes the mean value of the phase duration x, and σ2 is the variance. Note that, for the uniform distribution, the variability V is related to the standard deviation σ by the relation Embedded Image.

In the electronic supplementary material, figure S1 shows a typical time series obtained with the automaton model when considering lognormal distributions of the durations of the cell cycle phases. The mean durations of the phases are as in figure 5, and the standard deviation is equal to Embedded Image with V = 20%. In figure S1a, the inset illustrates the distribution of the duration of the S phase.

References

View Abstract