Dynamics in a delayed diffusive cell cycle model ∗

Abstract. In this paper, we construct a delayed diffusive model to explore the spatial dynamics of cell cycle in G2/M transition. We first obtain the local stability of the unique positive equilibrium for this model, which is irrelevant to the diffusion. Then, through investigating the delay-induced Hopf bifurcation in this model, we establish the existence of spatially homogeneous and inhomogeneous bifurcating periodic solutions. Applying the normal form and center manifold theorem of functional partial differential equations, we also determine the stability and direction of these bifurcating periodic solutions. Finally, numerical simulations are presented to validate our theoretical results.


Introduction
In organisms, the cell cycle, which processes DNA duplication and cell division, is the fundamental mechanism of cell proliferation.With plenty of biological experiments in the past few years, a large number of proteins, which are involved in cell cycle, have been revealed.These proteins form regulatory circuits that act like autonomous oscillators [15].The progression of cell cycle can be briefly divided into four phases: G1, S, G2 and M.During the G1 phase, the cell increases the organelle number and its size to prepare for the cell division.After that, the cell transits to S phase during which DNA is replicated.Following the S phase, the cell enters G2 phase and starts to synthesize some proteins, which are necessary for cell division.Finally, the cell passes the G2/M transition allowing the initiation of mitosis and divides into two daughter cells.The sequential occurrence of these four phases is driven by a family of cyclin/Cdks complexes [16,19].
Mathematical modeling is a powerful tool to investigate the dynamics of biological systems such as epidemiology [3,34], ecology [26,28,37,40] and molecular network [32,35].A number of mathematical models have been proposed to explore the regulatory mechanism of cell cycle in early embryos [10,18,21,31], yeast [6,7,27] and mammalian cells [1,9,22].Most of these models were described by nonlinear ordinary differential equations.These theoretical models revealed that the oscillation of a cell cycle system is driven by the delayed negative feedback loop provided by cdc2induced cyclin degradation.Furthermore, the bistable steady-state response originating from positive feedback loops was proved to enhance the oscillation robustness of cell cycle [20].Besides the comprehensive models mentioned above, analytical studies of cell cycle systems are also performed in some reduced models.Tecarro et al. used a 3-variable cell cycle model to describe G1/S transition [30].Using linear stability analysis, they obtained the stability conditions of the steady states in a simplified system.Moreover, they proved the occurrence of the periodic solutions emerging via supercritical Hopf bifurcations.
Recently, the spatiotemporal dynamics of cell cycle systems have been increasingly studied [5,33,38,39].Our previous work considered the diffusion of proteins in fission yeast and explored the effect of the spatial regulations on the initiation of mitosis [38].Vilela et al. constructed a reaction-diffusion-convection system interacting with a deterministic model describing microtubule dynamics and the cell size checkpoint of yeast [33].Ferell's group used partial differential equations to propose a theoretical model for revealing the role of bistability in the spatial activation of Cdks [5].Furthermore, the prediction is validated by their experimental data.However, most previous studies concentrated on the numerical simulations, and we haven't paid enough attention to mathematical analysis (the results on the stability and bifurcations) of the spatial effect on the dynamics of cell cycle systems.
The objective of this study is to analyze the stability of the positive equilibrium, the existence and stability of Hopf bifurcation induced by the delay in a cell cycle model involving spatial diffusion.The rest of the paper is organized as follows.In Section 2, we present a delayed diffusive cell cycle model for G2/M transition, which is system (3).In Section 3, the local stability of the positive equilibrium is obtained, and the existence conditions of homogeneous and inhomogeneous Hopf bifurcation periodic solutions are derived by analyzing the characteristic equations.In Section 4, the direction of Hopf bifurcation and the stability of the bifurcating periodic orbits are determined by using the normal form and center manifold theorem of functional partial differential equations.In Section 5, some numerical simulations are given to illustrate the theoretical results.Finally, we summarize our conclusion in Section 6.

The construction of a mathematical model
In this section, our aim is to establish a delayed diffusive model with Neumann boundary condition for depicting the spatiotemporal dynamics of the cell cycle system in G2/M transition.
Figure 1.Scheme for the mathematical model of G2/M transition.There is a delayed negative feedback loop in which MPF(u) promotes the activity of APC(v).In return, APC can reduce the concentration of MPF by increasing its degradation rate.In addition, MPF positively regulates itself via wee1 and cdc25.Thus, positive feedback loops are also involved in the model.
Figure 1 exhibits the core regulatory network in G2/M transition.MPF (mitosis promoting factor) is a kind of protein complex which promotes the initiation of mitosis.Cells can enter mitosis only when the concentration of MPF reaches a certain level.In this mathematical model, we assume that MPF is synthesized at a constant rate k s1 .Let u(t) represent the concentration of MPF.MPF can increase the production of APC (anaphasepromoting complex) represented by v(t).The basic synthesis rate of v(t) is denoted by a constant k s2 .In fact, the promotion of APC is not instantaneous, but mediated by some time lags required for intermediate biological processes.Thus a increasing function of MPF with time delay, g(u(t − τ )), is employed to denote the positive regulation of MPF to APC.In return, APC combines with MPF at a constant rate of k 3 and promotes MPF degradation, which closes a negative feedback loop [17].The constants k d1 and k d2 represent the basic degradation rate of MPF and APC, respectively.A pair of positive feedback loops are also involved in the regulation of G2/M transition (MPF activates its activator Cdc25 and inactivates its inhibitor Wee1) [8].In order to simplify the model, we do not show them in the figure and just use a increasing function, h(u), to represent the positive auto-regulation of MPF.Then, we can obtain the following delayed ordinary differential equations: where k s1 , k s2 , k d1 , k d2 , k 3 are positive constants.
Recent studies show that the diffusion of some proteins may play an important role in the dynamics of cell cycle [33].Experimentally, cell extract is usually loaded into a teflon tube to observe the diffusion dynamics.Thus, theoretical models investigating the spatiotemporal dynamics of cell cycle can be studied in one dimension.So we similarly introduce the one-dimensional spatial diffusion of MPF and APC into system (1).In addition, considering that MPF and APC in outside cells may not spread to the inside because of the cell walls, we use Neumann boundary condition for system (1).Therefore, a delayed diffusive cell cycle model with Neumann boundary condition is set up as follows: where d 1 > 0 and d 2 > 0 are the diffusion constants for MPF and APC, respectively.As an aid to perform analytical studies of system (2), we take h(u), g(u) as simple increasing functions h(u) = k 1 u and g(u) = k 2 u, where k 1 , k 2 are positive constants.For simplification of notations, we denote Thus, we obtain the following reduced model with Neumann boundary condition: where u(x, t) represents the concentration of MPF at time t and location x, v(x, t) represents the production of APC at time t and location x, a and c denote the basic synthesis rate of MPF and APC, respectively.f represents the basic degradation rate of APC.k 2 is the coefficient of the increasing function which denotes the positive regulation of MPF to APC. k 3 denotes the rate of APC combining with MPF.The difference of k 1 minus k d1 is denoted by b, and k 1 is the coefficient of the increasing function which represents the positive auto-regulation of MPF, k d1 represents the basic degradation rate of MPF.Here all the parameters a, c, f , k 2 , k 3 are positive, and b ∈ R. τ 0 is the time delay, l ∈ R + .
In the rest of the paper, regarding τ as a bifurcation parameter, we focus on the stability of the positive equilibrium and the Hopf bifurcation analysis of system (3).

Stability analysis of the model
In this section, by analyzing the associated characteristic equation, we investigate the stability of the unique positive equilibrium for system (3) and establish the existence of spatially homogeneous and inhomogeneous bifurcating periodic solutions depending on Hopf bifurcation.

Local stability of the model without delay
In this subsection, we obtain the stability of the positive equilibrium for system (3) without delay, which is independent of diffusion.Let τ = 0, system (3) is turned into Obviously, the equilibria of system (4) are determined by the following equations: Solving ( 5), we obtain two equilibria (u − , v − ) and (u 0 , v 0 ), where and It is not difficult to find that u − < 0, v − ∈ R and u 0 > 0, v 0 > 0, which show that only the unique positive equilibrium (u 0 , v 0 ) is in according with the reality.So we omit the non-positive equilibrium (u − , v − ) and just study the positive equilibrium (u 0 , v 0 ) throughout this paper.
For system (4) without diffusion, the Jacobian matrix at the positive equilibrium (u 0 , v 0 ) is It follows from ( 6) and (7) that Then, at the positive equilibrium (u 0 , v 0 ), the characteristic equation of system (4) is The roots of Eq. ( 9) are given by From (8), it is seen that the characteristic roots λ 1,2 < 0. So we have the following conclusion.
From the result of [23] one can see that the sum of the multiplicities of the roots of ( 12) in the open right half-plane changes only when a root appears on or crosses the imaginary axis.So we need to derive the conditions, which make the above cases happen.We first make the following hypothesis: and Then we have the following lemma.
In view of Lemma 1, concerning time delay τ , throughout the paper, we only consider the following case: Let λ n (τ ) = α n (τ ) + iβ n (τ ) be the root of ( 12) satisfying α n (τ j n ) = 0 and β n (τ j n ) = ω n when τ is close to τ j n .Then we obtain the following transversality condition.Then by ( 14) and ( 17), we have which means that https://www.mii.vu.lt/NADenote τ 0 * = min 0 k N0 {τ 0 k }.Based on the above analysis and the qualitative theory of partial functional differential equations [36], we have the following results on the stability and Hopf bifurcation.

Direction and stability of Hopf bifurcation
In this section, we discuss the direction and stability of Hopf bifurcation by using center manifold theorem and normal form theorem of partial functional differential equations [12,36].Similar approach has also been used in [14,28,40].
We first transform the positive equilibrium to the origin by the variable substitution ũ(x, t) = u(x, τ t) − u 0 and ṽ(x, t) = v(x, τ t) − v 0 .For simplification, we drop the tilde.Then for x ∈ (0, l) and t > 0, system (3) can be rewritten as follows: Fixing j ∈ N, for 0 n N 0 , we denote τ = τ j n .Let where L α (•) : C → X and G : C × R → X are given, respectively, by Consider the linear equation From Lemma 1 we can conclude that ±iω n τ are characteristic values of system (23), and the linear functional differential equation is as follows: By Riesz representation theorem, there exists a 2×2 matrix function η(σ, α), −1 σ 0, such that where and Let A(τ ) denote the infinitesimal generator of semigroup induced by the solutions of Eq. (24).
For a complete description of g 21 , we need to compute the expressions of W 11 (θ) and with then by using the chain rule Ẇ = W z ż + W z ż, we have Noticing that for θ ∈ [−1, 0), and when n = 0, we have In addition, for θ = 0, it follows from (28), (35) that where Then by (36) and the definition of A τ , we have for θ ∈ [−1, 0) where On the other hand, for θ = 0, it is concluded from the definition of A τ , ( 23)-( 25), ( 28) and ( 36) that Substituting ( 37) into (38), together with which implies that where n = 0, 1, 2, . . . .Furthermore, similar to the procedure of computing of W 20 , the expression of W 11 is as follows: where

Conclusion and discussion
A growing number of mathematical models have shown that the spatial diffusion may play a important role in the dynamics of various ecosystems [25], epidemic spread [2,4] and some genetic regulatory network.In this paper, we present a delayed diffusive model to investigate the spatiotemporal dynamics of the cell cycle system.Our theoretical analysis indicates that the diffusion may not play a crucial role in the generation of periodic oscillation.However, the time delay is more important, which determines the stability of the equilibrium: when the time delay is less than the critical value τ 0 * , the equilibrium is locally asymptotically stable.As τ increases and crosses τ 0 * , the equilibrium (u 0 , v 0 ) loses its stability and Hopf bifurcation occurs.Furthermore, we obtain a series of critical time delays, where the bifurcating periodic solutions are spatially inhomogeneous.
Bistability is a common phenomena of cellular function, which is involved in many biological systems such as cell cycle progression [5][6][7], cell differentiation [11,29], as well as the development of cancer [13,24].Although we do not consider the bistability, which originates from the positive feedback loops, in our model, the results in [5] show that the bistable steady states may accelerate the spatial activation of Cdks in Xenopus cell cycle.Also, trigger waves can be observed experimentally.Therefore, we will explore the effect of bistability on the spatiotemporal dynamics of cell cycle systems in the future work.

1 , 2 ,
. . ., are the eigenvalues of ∆ in [0, l] with the Neumann boundary condition.It follows from (11) that the characteristic equations at (u 0 , v 0 ) are the following sequence of quadratic transcendental equations