Firstly, since the purpose of this model is to explore the dynamics of the cell cycle within a single-cell, we describe the ‘state’ of the cell based on the relative amounts of the molecules which govern this process; i.e., the total mass of cyclins at time t. For example, the S-cyclin state of the cell describes the relative amounts of S-cyclin in the cell, whilst the overall state of the cell is determined by the relative amounts of S- and M-cyclins or the lack thereof. This state then allows one to determine both the phase of the cell cycle, within which the cell finds itself, and also the trajectory of the cell in time.

Now, let us consider the S-phase of the cell cycle. Allow that all dynamics of this system occur within a time interval, I:=[0,T], and an internal cellular spatial domain, D⊆ℝd, of dimension d with independent variables t∈I and ξ∈D (Table 1). Since we are assuming that only S- and M-phase dynamics are relevant for our system, the expression of S-cyclin will necessarily be tied to the activation of the M-phase in the previous cell cycle, such that the S-cyclin state should not increase whenever there are significant numbers of M-cyclins still present within the cell.

In order to account for the change in the bio-chemistry of the cell, without introducing additional states, we then allow for the existence of a “cellular memory”, whose length is determined by the temporal parameter τ. Biologically, cellular memory can accumulate through the biochemical transformation of a protein's surroundings (for instance, through phosphorylation of target proteins), so that the biochemical landscape of the cell will retain some residual differences from its base state, even after the disappearance of the protein [40, 41]. In this case, and in order to reduce the number of independent variables in the model, we account for this by allowing a delay-intergral in the system of integro-differential equations. In practice, τ serves to determine the time-scale over which the cell may account for its past states while, in theory, this τ represents the time-scale over which the state of a cell will determine its future state, through influencing biochemistry in the present. Thus, we assume that if the historically weighted M-cyclin state, giving rise to current protein levels, is below some critical threshold, θ_m, then it is assumed that S-cyclin expression will occur, and not in any other case.

Moreover, to account for the fact that cyclins will be produced nearby the nucleus and must diffuse throughout the cell to control whole-cell cell-cycle functioning, we multiply this positive term by appropriate spatial biasing kernel p:D→ℝ+. The purpose of the spatial biasing kernel is to distribute the translation or production of cyclin proteins in the region near to the nucleus and, therefore, we choose the form

Where κ_p is a decay rate from ξ = 0. Therefore, assuming that the S-cyclin state is normalized and asymptotic to its maximum quantity, the term corresponding to the expression of S-cyclins (i.e., “source of S-cyclins”) within the cell is written as

Where the positivity operator is given by [·]⁺ : = max(0, ·) and k_s > 0 is the cyclin production rate constant.

In simple terms, when m(t, ξ) is sufficiently high across space and some period of time and s(t, ξ) < 1 for p(ξ) > 0, the source term should be positive and S-cyclins are produced. If either of these conditions are not met [i.e., either m(t, ξ) was not sufficiently high for long enough or s(t, ξ) > 1] then the source term will be 0. This acts as a switch to engage production of cyclins and begin the given phase of the cell cycle, just as in the biological case.

Further, the degradation of S-cyclins and subsequent decrease in amount appears, also, to depend upon the M-cyclin state, with S-cyclin being systematically degraded as M-cyclins approach their maximal cellular quantity. Therefore, assume that degradation occurs with a rate μ_s, we have that the total S-cyclin decay is given by

where

is the usual characteristic function.

On the other hand, exploring now the balance between the source and decay of the M-cyclin population, the relationships between the involved positive and negative terms (of identical form to the ones giving the source and decay for S-cyclins) is here reversed, as M-cyclins are expressed only when historic S-cyclin population levels are above a certain threshold and are degraded only when current S-cyclin states approach their minimum value, 0. Furthermore, both for S- and M- cyclin populations, given that decay toward minimal cyclin quantity is regarded as an exponential asymptotic process, explore this by introducing a parameter ϵ which represents the degree to which levels of cyclin are expected to approach their respective maxima or minima.

Finally, in the presence of the associated balance of source and decay, in order to permeate the entire cell, per unit time, both cyclins populations are assumed to exercise a diffusive spatial transport locally throughout the cyctoplasm, with diffusion coefficients D_s and D_m for S- and M-cyclins, respectively. Hence, the coupled dynamics of S- and M-cyclin populations is therefore given by

Given that these dynamics are assumed to take place within the cyclin-impermeable lipid membrane of the cell and, as such, we impose no flux boundary condition to complete the model.

Beginning with initial conditions given by s(0, ξ) = m(0, ξ) = 0 and using a finite differences numerical simulation scheme, alongside central difference calculations for gradient approximation and mid-point integration calculation for integrals, we obtain both spatially and temporally resolved results. As such, we present these firstly as a spatially integrated cyclin state profile over time (Figure 1) and secondly as a time series distributed in the internal cell space, ξ (Figures 2–4). The default parameter set used in this section is given in Table 2.

Observing the temporally resolved data, in the first instance, we find an oscillatory dynamics between the S- and M-cyclin states (Figure 1B). We also observe that the time period over which S-cyclin maintains its near-maximal quantity, |1 − s| < ϵ, is significantly longer than that observed for the M-cyclin, in spite of the controlling variables being equal in value: θ_m = θ_s. Also, the relative exclusivity of these two cyclin states means that we can begin to delineate the S- and M-phases, respectively, as those time periods across which either the S- or M- cyclins are at their maximal concentrations, yielding a period of Δt = 2,770.

Let us turn, now, to the spatially resolved system, wherein the S- and M-cyclins are viewed as a time series, as they diffuse through the spatial cell interior, ξ∈D. Beginning at t = 10, we observe the creation of S-cyclins at the origin of the system, ξ = 0, consistently with the cyclin translation (or production) function, p(ξ), followed by the asymptotic homogenisation of the spatial domain through diffusive processes for t ∈ (0,500] (Figure 2). At t ≈ 1,000, the historic activity of S-cyclins becomes sufficient to manifest the process of M-cyclin creation which, again, begins from the origin and homogenizes across the domain as time progresses (Figure 2), such that S- and M-cyclins coexist for some period of time.

At t ≈ 1,310, the population of S-cyclins decreases ubiquitously across the domain as ∫Dm(t,ξ)dξ′+ϵ>1 (Figure 3). It may also be readily observed that this process is exponential and asymptotic with s(t, ξ) → 0. Likewise, at t = 1,770, the M-cyclin population begins to decrease ubiquitously in space and exponentially in time, such that m(t, ξ) → 0.

Comparing this spatial process (Figures 2, 3) with its corresponding temporal process (Figure 1B), we find a close correspondence where the diffusive process accounts for the delay between process initiation and concentration maximization in cyclin states. Moreover, observing the temporal sequence, we find that the S-cyclin state remains significantly greater than 0, whilst the M-cyclin state appears to be at least infinitesimally close to 0. (In fact, since the processes of protein production and protein destruction are both asymptotic, the M-cyclin states may be no less than infinitesimally greater than 0 for times, t′, at which m(t, ξ) > 0 with t>t′,ξ∈D). This is due to the assumption that S-cyclins are degraded only for times in which M-cyclins are near-maximal, whilst M-cyclins are degraded where S-cyclin states are near-0, inducing a fundamental asymmetry between these quantities; the time period over which each are degraded is significantly different, under these assumptions and with the given parameter values.

Finally, as the M-cyclin state remains sufficiently low for a sufficient period of time, t ≈ 2, 770, one observes the S-cyclin state begin to rise (Figure 4), similarly as in the above case (Figure 2, t = 10). Likewise, at t ≈ 3, 740 the M-cyclin state again begins to rise and diffuse across the cell to complete the periodic oscillation (Figure 4), correspondingly with the temporal cyclin state data (Figure 1B).

Since the derivation of the G-phase dynamics is analogous to that of the S-phase, we begin here simply by stating the fundamental dynamical equation in the G-phase:

As such, the full derivation of the dynamical equation for the G-phase follows that given in Section 3.1.1.

Since the derivation of the M-phase dynamics is analogous to that of the S-phase, we begin here simply by stating the fundamental dynamical equation in the M-phase:

As such, the full derivation of the dynamical equation for the M-phase follows that given in Section 3.1.1.

Therefore, with all of the individual elements of the model introduced in Equations (10)–(26), the full system of equations for the 3 phasically-resolved subpopulations are given by:

Now that the system, in its entirety, has been derived and discussed, it is useful to consider how various parameters will affect the behavior, dynamics, and outcomes of the system. Firstly, an increase in synthesis rates (χ_s, χ_g, and χ_m) will cause a concurrent increase in the dynamic transition rate of cells to the subsequent phase of the cell cycle, as a result of their accumulation of cyclins and favorable positioning with respect to the relevant ψ. Likewise, an increase in random cyclin fluctuations (D_s, D_g, and D_m) will both cause some cells to defer transition to the subsequent phase and others to transition more rapidly, as a result of the widening distribution in cyclin concentrations.

Now consider the more complicated dynamics of the G2-phase: Firstly, the p53 concentration, p, will increase only with lower levels of p53 and mdm2, μ, since its factor (1-μωμ-p) is near-zero for high levels of mdm2. Likewise, the mdm2 concentrations will increase only for low levels of DNA-damage and mdm2, since (1-μωz-μ) is near-zero for high levels of DNA-damage. The production of mdm2, however, is also contingent upon the presence of p53. Thus, for high levels of DNA damage, p53 will be produced in the absence of mdm2 until the DNA-damage has been sufficiently repaired, with rate χ_z, to allow for increases in mdm2 concentration. At such times, p53 production will reverse to a degradative dynamic to allow for p53-dependent transition, with θ(p). Thus, increasing the rates χ_p or χ_μ will significantly effect the extent to which this checkpoint operates. For instance, with low χ_p, one would expect high rates of transition from DNA-damaged cells, due to the lack of p53 to suppress transition. Increasing χ_μ will likewise admit greater numbers of cells to the subsequent phase due to the increased rate of suppression of p53 in the G2-phase. Consider, also, that an increase in either of the diffusion parameters in the DNA-damage space, D_z and D_r, will therefore cause some subset of cells to remain within the G2-phase for longer, on account of this mdm2-p53 checkpoint.

These inter-phase dynamics yield interesting results for which the trend may be approximated, though the outcome of simultaneously altering multiple parameters may be difficult to predict, a priori. We can also understand general principles of the system, such as that having a greater ϕ_g, at the expense of ϕ_s (since ϕ_sϕ_g = 2) will allow for a greater efficiency in the system, since biosynthesis is not unnecessarily spent on cells with extensive DNA damage.

In the first instance, we began by simulating the cell cycle by generating solutions to Equation (27), using the given parameters (Table 3), and extracted population level metrics during simulation in MatLab. As expected, cells lost from the 1st cell cycle are incorporated into the second cell cycle, while each cell cycle population exhibits a sigmoidal growth phase before and apparently sub-exponential decay phase (Figure 5A). It is also worth noticing that twice as many cells enter the 2nd cell cycle as have left the 1st so that the population is doubled through a successful mitosis.

Since the processes of the cell cycle are probabilistic, rather than deterministic, and rely on several checkpoints in order to complete a full cycle, the 1st cell cycle still retains a significant percentage of its initial population, even at t = 150. This is despite a speedy and almost complete transition of cells from the S-phase with the time interval t ∈ [0, 100]. In developing organisms, such as Drospohila melanogaster, it is therefore likely that cellular biochemistry is more tightly regulated, in order to achieve the tightly regulated progression of the cycle across the population and the synchronization of mitosis across the embryo during early cycles.

As for the phases of the cell cycle, the entirety of the population originates in the S-phase of the cell cycle, by design (Figure 5B). Although this constitutes a departure from the biological reality of the average system, it allows for a more systematic observation of dynamics and divergence across the system, as a whole. As cells leave the S-phase of the cell cycle, they enter the G2-phase of the cycle and quickly begin to progress through to the M-phase. It is in this G2-phase of the cell cycle that the individual cells, amongst the population can be seen to diverge in their behaviors. Although there exists a significant transfer of cells between the S- and G2-phases, this progresses is arrested in the G2-phase due to the dual requirements of the cell to achieve low levels of DNA-damage and high levels of G-cyclin. This leads to an early onset but slow rate of transfer between the G2- and M-phases of the cell cycle. Cells which do enter the M-phase, however, appear to progress through mitosis successfully and reach their next cell cycle.

Similarly as what is observed in the case of numbered cell cycle populations (Figure 5A), the population is seen to increase at around t = 30 so that the doubling time for the population is given approximately by 150 (Figure 5C). Interestingly, however, the increase in the population has a linear trend, where one might have expected an exponential trend for such a system. Looking back toward the populations within various phases of the cell cycle (Figure 5), it is possible that this linear trend in population growth is explained by the apparently linear trend in the growth of cells within the M-phase, owing to G2-phase checkpoint controls. It is difficult to know whether this relationship does indeed hold in nature or whether the rate of these molecular checkpoints allows for an exponential realization. In order to establish the mass-balance of the cell cycle decomposition approach, we have provided an overlay of the integral across the cell cycle phases, giving the total population, and the sum of the populations across the relevant cell cycles.

Of course, as is the case in the CDK cell-cycle model, the cell-cycle is characterized by the transition from the S- through the G- to the M-phase biochemistry, where each of these phases are themselves characterized by unique dynamics and cellular behaviors which ultimately result in the growth of the cell population. The progress of the studied cell population through the cell cycle is seen from the decompositions given in Figure 5A, whilst the dynamics of the solutions for the cell population are found in Figures 6–11.

It is important that one interpret these dynamics carefully, so as to avoid leaving with the mistaken impression that the true nature of the population's dynamics are given explicitly by these solutions to system (Equation 27). Instead, one must bear in mind that these solutions describe density distributions for the population, which must, in turn, be interpreted as a distribution of possible states in which one might expect to find the system. It is the sum of probabilistic expectation functions, across cells, so that each value is the multiple of a population-wide probability density function and the total number of cells at time t. Thus, when small densities are first recorded within the G2 phase, this implies only that there is some non-zero probability that a cell has (or cells have) transitioned into the G2 phase.

As in the case of demographic distributions across a given human population, for example, these density distributions describe the likely characteristics of the populations of a whole, whilst they may fail to describe the precise individual realizations of such a process. For instance, outlier events occur frequently and confound the rationale for such a demographic distribution and, yet, the distribution remains an accurate depiction of the large-scale populations dynamics. Similarly, we here depict the dynamics of the density distribution for entire population, such that fractional elements of cells may transition between phases of the cell cycle, truly representing only a fractional probability that a single whole-cell may have achieved transition across such a boundary. One must appreciate, therefore, that what may be true for the population may still yield surprising results on the individual scale.

Next, in order to investigate the effect of reducing the rate of cell population's intrinsic DNA repair pathway, we re-parameterise the system for χ_z = 0.02 and observe the change in the underlying structural dynamics in comparison to the above case (Section 4.2.2). During the early stages of simulation (t < 10), results remain qualitatively and quantitatively similar to previous cases and, as such, these results have been omitted for the sake of scientific parsimony.

As significant numbers of cells begin to accumulate within the G-phase of the cell cycle (t ∈ [10, 30]) one observes a significantly increased standard deviation in the DNA-damage spectrum of the cellular population (Figures 12, 13) when comparing to the equivalent population in the previous case (Figures 9, 10). This is to be expected as a result of reducing the DNA-damage repair rate, since χ_z controls the rate at which the cellular population exhibits asymptotic advection toward z = 0, by definition. On the other hand, the currently considered cellular population (Figures 12, 13, Cgp,μ) exhibits a similar rate of increase in automatic p53 expression across the early stages their G-phase occupation, compared with the previous case (Figures 12, 13, Cgpμ). Mdm2 expression, as a result of the reduction in DNA-damage repair rate, is, however, limited in this case and transforms the distribution from its previous quadratic form (Figures 12, 13, Cgpμ) to a far more shallow distribution, exhibiting a form of an order 2 > q > Q polynomial (Figures 12, 13, Cgpμ), where Q > 2. We must bear in mind that this is a distribution, such that cells exhibit multiple positions upon the (p, μ)-plane which lie between the curves bounded by the order 2 and order Q polynomials.

Simultaneously, we begin to observe an increase in the p53 expression levels (Figures 12, 13, Cgp,z), previously not exceeding p ≈ 0.5 (Figures 12, 13, Cgp,z), as the DNA-damage levels decrease at a slower rate. During later time points, t ∈ [70, 150], one observes a significant increase in the p53 expression levels in cells with increased residual DNA-damage (Figures 13, 14, Cgp,z) in comparison cells within the previously considered population (Figures 13, 14, Cgp,z), who typically underwent DNA-damage repair prior to the occurrence of this increase in p53 expression. Thus, the cellular population appears to mitigate this additional residual DNA-damage by restricting the passage of the cellular population between the G- and M-phases at the G-phase checkpoint.

Once again, the cellular distribution for the population with reduced DNA-damage repair rate exhibits a significantly differing profile in the (p, μ)-plane (Figure 13, Cgp,μ) when compared to the previous case (Figure 10, Cgp,μ). In the current case, the distribution is not, in fact, symmetric whilst the previous case shows a reflective symmetric distribution in μ, whilst the additional qualitative comparison on these 2 distributions is left to the reader. Similarly as with the previous case (Figure 11, Cgp,z), however, we observe that the cell distribution which migrates negatively in p has considerably lower mean DNA-damage levels when compared with the population as a whole (Figure 14, Cgp,z).

Likewise, across the whole time sequence (Figures 12–14, Cgg,z), the increase in G-cyclin remains consistent with the previous case (Figures 12–14) whilst the DNA-damage remains elevated, with respect to the standard deviation of the distribution. Notice, further, that the modal DNA-damage is equal in both case, at z = 0, since the DNA-damage mechanism is purely diffusive (such that there exists no mechanism by which this modal value could change). Likewise, even the selective transfer mechanism between the G- and M-phases should ensure that the modal DNA-damage rate within the transferred population remains z = 0. Therefore, to discuss dynamics in the DNA-damage dimension in the case where the initial condition has a modal density at z = 0 is always to speak of the change in the populations distribution through z and not its position in z, per se. This means that although individual cells may change their position with respect to DNA-damage, the most likely situation for a cell is remain without DNA-damage throughout the cell cycle.