The important role relevance between modelling phase transition and complex dynamics of the cell cycle is found beyond the knowledge of biology background yet. The reason lies in that, on one respect, mathematical models of cell cycle can contribute to the understanding of its mechanism biologically; On another respect    , the dynamical analysis mathematically is developed as a ubiquitous technique to predict system evolutionary behavior and the results provide the testable suggestions. For example, the robust stability of the state of the “checkpoint” in G0/G1 phase .
In mathematical description of biological models, some key components and their mutant relationship with circumstance (inner-cell or inter-cells) are considered. People always translate their necessary biology knowledge into some differential equations to describe the growth stage of during cell cycle which include different phases   . As is well known, G0/G1 “gaps” joint the S phase and M phase to form the cell cycle. Hence, we introduce time delay to model the necessary time experienced to transpass “gaps”. Time delay is incorporated into the phosphate groups which while binding to target proteins to active protein phosphorylation process. In addition, to drive the downstream events to generate gene protein production, enzymes reaction formed by binding itself to Cdc2 monomer which is assumed no degradation in cyclin synthesis  . Hence, Cdk2 dimer dynamics is dominated by the well-known Michaelis-Menten rate laws. The Michaelis-Menten rate law is described as
In cyclic synthesis, the rate of association with Cdk2 monomers depends on the probability of a collosion between two species. The enzyme activation/suppression leads to target protein production generated to form cyclin. Therefore, the protein phosphorylation and dephosphorylation activity dramatically altered the target protein which regulates kinase of cell growth.
Time delay describes the present state or state feedback dependent on the past history which becomes differential equations into DDEs of infinite dimensional  . Based on the above biology background, we set up delay differential equations (DDEs) to express protein regulation in the cell cycle, and write G1 mitosis phase into a useful mathematical expression.
Besides, one mother cell can be divided into two daughter cells during mitosis stage; the dynamics in G1/M phase can also be expressed as
wherein x denotes enzyme concentration, and y represents gene protein concentration, is time delay with its valuable estimation between 2.8 h and 7.8 h. The associated function is written as
State feedback control is to introduce concentration difference into the model which reflects the time delay effects on dynamical evolution behavior. Herein , with denotes time delay in “gap” junction as G0 phase.
System (1.2) is DDEs with multiple time delays which induce complex dynamical phenomena underlying Hopf bifurcation, such as the coexistence of four periodical solutions. It is observed that Neimake-Sacker bifurcation phenomena arise as varying time delay, and the quasi-periodic solution is also observed in coexistence regime of periodical solutions. By applying DDE-Biftool technique  , the continuation of periodical solution bifurcating from Hopf bifurcation is carried out by varying time delay. The exhibited limit cycle dies out in the vicinity of unstable region partitioned by Hopf lines from stable region. The hysterisis phenomena of steady state are observed and “Z-shaped curve is pictured since the collision of stable equilibrium and unstable equilibrium at points with fold singularity.
The saddle-node bifurcation of limit cycle occurs at some critical values of time delay . The stable periodical solutions are continued though the instability oscillation of bifurcating periodical solutions birth from both sub/super-critical Hopf point. The hysteresis phenomena of limit cycle are also observed beyond Hopf bifurcation. As is well-known, based on the strong continuous semigroup theory, the linear version of DDEs is described by the related form of ODE operator. Hopf bifurcation occurs as the infinitesimal generator has a pair of simple imaginary roots. A possible approach to investigate this bifurcation problem for DDEs involves the computation of normal form on center manifold. The preceeding work of normal form approach is famous with Faria and Magalhaes’s work, see also      for reference. The coefficients of the normal form are computed by the Lyapunov-Schmidt dimensional reduction method, and the formula for calculating coefficients is also expressed in Section 4   , and the formula is applied in Section 3 directly to compute the bifurcation direction and stability of bifurcating periodical solutions of Hopf bifurcation.
The whole paper is organized as follows: In Section 2, the bistable dynamics of steady states are analyzed. In Section 3, the imaginary roots of the related characteristic equation of the linear version of system (1.2) are analyzed to derive the critical values of Hopf bifurcation. With the introduction of a new time delay and define three different angle variables, the threshold of Hopf bifurcation lies in the intersection of two parameterized curves. The bifurcating periodical solution is continued as varying time delay continuously. In Section 4, the normal form near Hopf point is computed which further discovers that the bifurcating periodical solution is unstable, and this further verifies the numerical simulation results obtained in Section 3.
2. Bistability Phenomena in G1 Phase of Cell Growth
The bistability phenomena is observed as varying free parameter pairs to track fold bifurcation curves of the equilibrium solution. The cusp degenerate singularity is also tracked as the codimension-2 bifurcation point whilst fold curves on parameter plane immerge, as shown in Figure 1(a) and Figure 1(b). We define the curve
Figure 1. The equilibrium solutions of system (1.1). (a) The bistability property of system (1.1) is observed since two LP bifurcation points exist; (b) The codimension-1 bifurcation as fold bifurcation curves are plotted, and the immergency of the two-fold curves leads to cusp sigularity; (c) The unstable region of equilibrium solution of system (1.1) is partitioned from the stable region( inner region and outer region of Γ curve).
As shown in Figure 1(c), the region enclosed by the curve is unstable region lying on three-dimensional space, and the curve partition the outer stable region from the inner unstable region. It is observed that one stable equilibrium solution collides with the unstable one as the parameter pair is attached to the vicinity boundary of fold bifurcation curves. The cusp bifurcation point CP is denoted as , system (1.1) exhibits bistability property if choosing free parameter less than .
3. Hopf Bifurcation in G1/M Phase
Hopf bifurcation triggers the phenomena of small amplitude periodic oscillation. Hopf bifurcation happens as one pair of characteristic roots cross the imaginary axis transversally as varying free parameter and time delay continuously. To study Hopf bifurcation, people used to linearize system near the positive equilibrium solution and further to analyze the roots with zero real part of the related characteristic equation.
In system (1.2), we assume the positive equilibrium solution is which satisfies
We plot x-nullcline and y-nullcline into x-y plane, the intersection of two nullclines is exactly the value level to describe the magnitude of the equilibrium solution, as shown in Figure 2. The linearized equation at is given as
Therefore, the characteristic equation of the linearized version (3.2) is written as
Figure 2. Hopf bifurcation of the positive equilibrium solution. (a) The continuation of equilibrium solution as varying free parameter. Hopf bifurcation occurs at or , other parameters are fixed as , , , , , , , , , ; (b) The curve of angle and with respect to the parameter is plotted with . The intersection point is associated with Hopf bifurcation; (c) The curve of and respectively with respect to is plotted with ; (d) The curve of and respectively with respect to is plotted with .
Furthermore, one obtains that
where is the free parameter.
Suppose , substitute it into Equation (3.5) and separate the real part from the imaginary part to get
with represents the real part and the imaginary part of the polynomial respectively.
Solving from Equations (3.6) to get
Since in DDEs, some delays are specified which manifest very difficult dynamical properties as its corresponding characteristic equation is a transcendental equation. With the introduction of a new time delay, which is not really represented in DDEs, however, though undefined originally, a new time delay can be applied to solve unknown time delays. Accordingly, by introducing new time delay and the associated angle variables , we define the new mapping
with , then rewrite the polynomial as
Define the new function
By Equation (3.10), one defines curve L by
By Equation (3.11), we also get the description of curve S
In geometry, Equation (3.11) and Equation (3.12) exhibit the infinite intersection or no intersection of curves L and S, and we list Hopf bifurcation condition as and for some k and l. With the assumption that the intersection of two curves is expressed as which satisfies , we compute the time delay as the critical value of Hopf bifurcation, . Note that is represented by the introduction delay s. To calculate the Hopf bifurcation line, Equation (3.11) builts the inverse function with lying inside some subinterval of since we have defined the relationship between and by mapping (3.8), and Hopf bifurcation occurs at
Following we compute the transversal condition of Hopf bifurcation, differentiate the rightside of Equation (3.4) with respect to , one obtains
then by solving from Equation (3.14), we get
differentiate the rightside of Equation (3.4) with respect to to further get
Noticed that and , we get , then substitute it into Equation (3.16) to get
then by Equation (3.17) and Equation (3.15), one has
herein denotes the conjugate part of , and expresses the conjugate part of . Noticed the expression of the characteristic polynomial (3.4), one has
The stability property of the equilibrium solution changes into unstable state as ascending free parameter and stability switching phenomena occurs, as shown in Figure 2(a). Respectively, we choose to further continue equilibrium solution by varying time delay and the stability switching phenomena are observed, as shown in Figure 3(a) and Figure 3(c). Hence the intersection of curve L and
Figure 3. Hopf bifurcation as varying time delay and Hopf line on parameter plane. (a) The stability switching phenomena for the equilibrium solution with chosen ; (b) Hopf line on parameter plane with ; (c) The stability switching for the equilibrium solution as ; (d) Hopf line on parameter plane with .
S satisfies , as shown in Figure 2(b). Hopf bifurcation occurs which specifies angle parameter to derive time delay . The curve S and L are drawn in plane. By equality , the bifurcation threshold of Hopf bifurcation is approximated at the intersection of two curves, as shown in Figure 2(c) and Figure 2(d). As the above description, the stability switching phenomena arise as varying rime delay , as shown in Figure 3. In Figure 3(b) and Figure 3(d), Hopf lines are drawn on parameter plane. Using DDE-Biftool software, the coexistence of four periodical solutions is shown in Figures 4-6 in Section 4. To compute the stability of the bifurcating periodical solutions with small amplitude, the coefficients in normal form are computed with the aid of dimensional reduction method which is given in Section 5. We have computed that the bifurcating periodical solutions are unstable near Hopf point and Hopf bifurcation is subcritical.
Numerical Simulation of Limit Cycle Continuation
Varying free parameter , Hopf bifurcation occurs at and , and periodical oscillation phenomena are observed which arise at Hopf points. By applying DDE-Biftool software, the bifurcating periodical solutions are calculated with high technique and which further simulate the periodical
Figure 4. The coexistence phenomena of four periodical solutions. (a) near ; (b) near ; (c) near ; (d) near .
Figure 5. The phase portraits and time series solutions of four periodical solutions as denoted by 1, 2, 3, 4 in Figure 4. (a) The phase portraits of four periodical solutions as numbered by 1, 2, 3, 4 in Figure 4(a); (b) The corresponding time series solutions of four periodical solutions; (c) The phase portraits of four periodical solutions as numbered in Figure 4(b); (d) The corresponding time series solutions;
solutions continuously as varying free parameter or time delay respectively. As shown in Figure 4(a) and Figure 4(b), the coexistence phenomena of bifurcating periodical solution are observed as varying free parameter near subcritical Hopf point or continued by supercritical Hopf point respectively. The red line represents the unstable periodical solutions while the blue line denotes the stable periodical solutions. The coexistence phenomena of multi-periodical solutions as varying time delay are also shown in Figure 4(c) and Figure 4(d). By choosing four points in Figures 4(a)-(d), the corresponding periodical solutions are also shown in Figure 5 and Figure 6, wherein the unstable periodical solution is drawn by red color or purple color, whilst the stable periodical solution is drawn by blue color or green color. In Figure 5 and Figure 6, the phase portraits of four periodical solutions and the corresponding time series solutions are pictured colorfully.
Hopf bifurcation brings forth the stability switching phenomena of the equilibrium solution of system (1.2). On one respect, Hopf bifurcation partition the stable regime from the unstable regime hence the stable equilibrium solution
Figure 6. The phase portraits and time series solutions of four periodical solutions as denoted by 1, 2, 3, 4 in Figure 4. (a) The phase portraits of four periodical solutions as numbered in Figure 4(c); (b) The corresponding time series solutions; (c) The phase portraits of four periodical solutions as numbered in Figure 4(d); (d) The corresponding time series solutions;
changes to be unstable when crossing over the margin of the stable region then switches back to be stable state again as further crossing over the margin of the unstable region. In addition, both the period-2 oscillating solution and the quasi-periodical solution manifests in system (1.2) as ascending time delay, as shown in Figures 7(a)-(d) and Figures 8(a)-(d). By Poincare section analysis, the switching phenomena between the period-2 solution and the quasi-periodical solution are also observed by numerical simulation technique, as shown in Figure 9(a) and Figure 9(b). It is remarked that quasi-periodical solution arises corresponding to the stability switching phenomena of the equilibrium solution, and the red circle denotes the corresponding Poincare section in Figure 7(d) and Figure 8(d).
4. Center Manifold Analytical Method
It is verified that Hopf bifurcation occurs at parameter pairs as shown in Figure 3(b) and Figure 3(d), since the property of the characteristic roots with zero real part. As is well known, by applying center manifold analytical technique, people always obtain the reduction dimensional system near Hopf point by Schmidt-Lyapunov scheme. However, it is believed that parameter perturbation always leads to a direct method to acquire the classifying technique by varying small coefficients in formal norm continuously. Based on parameter perturbation method, we compute the direction of Hopf bifurcation singularity by normal form computation and further to analyze the property of stability of the bifurcating periodical solutions.
With the introduction of , set , and scale , , then system (1.2) is transformed into its 3rd Taylor expansion truncation
Figure 7. The time series solution and phase portraits of system (1.2). (a) The time series solution of the period-2 solution with parameter ; (b) The corresponding phase portraits of the period-2 solution; (c) The time series solution of the quasi-periodical solutions with ; (d) The phase portraits of the quasi-periodical solutions.
Figure 8. The time series solution and phase portraits of system (1.2). (a) The time series solution of the period-2 solution with parameter ; (b) The corresponding phase portraits of the period-2 solution; (c) The time series solution of the quasi-periodical solutions with ; (d) The phase portraits of the quasi-periodical solutions.
Figure 9. The routes of switching phenomena between the period-2 solution and the quasi-periodical solutions. (a) The Poincare section with ; (b) The Poincare section with .
with initial value definition
wherein . According to the property of DDEs, solutions of Equation (4.1) are continuous on Banach space . With few discontinuity jumps, the solution operator is differentiable on the extended phase space . With the assumption , and denote for any , then Equation (4.1) is written as
Based on the fundamental theory of DDEs, there is a bounded variational function to express the linear version of Equation (4.2) on the extended phase space. In addition, by Reize representation theorem, the corresponding 2 × 2 bounded matrix function is listed below to write the linear version into its integral operator form, that is, for any ,
and the nonlinear part is written as
Consider the linear operator (4.3), it’s an infinitesimal generator of the strong continuous semigroup in phase space BC, and for any , we define new linear operator and its adjoint operator wherein , that is,
For , define the bilinear form
Based on the above analysis, Equation (4.2) can be written as its operator differential equations, that is
Considering the linear version of differential operator (4.8), Hopf bifurcation occurs at , and the associated characteristic roots set is denoted as as verrfied in Section 2. With the assumption of other eigenvalues with negative real parts, the phase space can be suspended on its restricted center manifold. Hence the eigenspace is decomposed into the eigenspace P corresponding to roots set and its complementary subspace Q. We suppose eigenspace P is spanned by where and its conjugation vector are respectively the eigenvector of and . We represent the eigenbasis of associated with linear operator and its adjoint operator respectively, then
By the definition of Equation (4.8), we have .
With a possible discontinuity jump at , we also define the mapping on the phase space as the following
for any , we can write it as , define the projection mapping , then we have , therefore, we can write as the expression with . Further computation express the differentiate relationship of axis variable given as
Set , then Equation (4.12) is written into its Taylor expansion with 3rd truncation to be expressed as
Define the operator on the homogeneous space as
Then the normal form on the center manifold can be written as
Hence by (4.15), one obtains the formal norm near Hopf point
The bases of projection are easily to compute as
The bases of are expressed as
Hence, we get
we eliminate the term . Set
then by Equation (4.13) and Equation (4.14), for , we have
The initial condition of Equation (4.22) is also computed as
Based on the initial condition (4.23), one can derive the coefficients by the integral Equation (4.22) with respect to . Therefore, we get the coefficient in formal norm (4.18) to be listed as
Cell grows in size by double in mass then divide to form into two daughter cells. In such a way, cells segregate to form more descendant cells to make lifetime perpetual both in plants and animals. The mutation relationship between different reaction components during cells growth and signal transaction between phase transitions is ubiquitous in regulating cells growth and division. The introduction of time delay is a crucial factor in cells cycle growth model since signal transaction is experienced in G1/G2 phase. The cells growth model was set forth by the enzyme reaction formed by binding itself to Cdk2 monomer to drive the downstream events to generate gene protein production. In addition, the protein phosphorylation and dephosphorylation activity dramatically altered the target protein which regulates kinase of cell growth. We focused on the stability analysis and bifurcation phenomena in cells growth system. Hopf bifurcation occurs as varying free parameter and time delay continuously. The co-existence periodical solutions were observed as varying parameter continuously near Hopf bifurcation point. The manifest stability switching phenomena also brings forth the quasi-periodical solution appearing in system. The complex dynamics as switching routes between period-2 solution and quasi-periodical solution were explored. We applied center manifold analytical scheme combined with the dimensional reduction method to analyze the normal form near Hopf bifurcation, which further verifies that Hopf bifurcation is sub-critical and the bifurcating periodical solution is unstable.
 Novak, B. and Tyson, J. (1993) Numerical Analysis of a Comprehensive Model of M-Phase Control in Xenopus Oocyte Extracts and Intact Enbroys. Journal of Cell Science, 106, 1153-1168. https://doi.org/10.1242/jcs.106.4.1153
 Angeli, D., Ferrell, J.E. and Sontag, E.D. (2004) Detection of Multistability, Bifurcations, and Hysterisis in a Large Class of Biological Positive-Feedback Systems. Proceedings of the National Academy of Sciences of the United States of America, 101, 1822-1827.
 Tyson, J. (1991) Modelling the Cell Division Cycle: Cdc2 and Cyclin Interactions. Proceedings of the National Academy of Sciences of the United States of America, 88, 7328-7332.
 Goldbeter, A. (1991)A Minimal Cascade Model for the Mitotic Oscillator Invovling Cyclin and Cdc2 Kinase. Proceedings of the National Academy of Sciences of the United States of America, 88, 9107-9111. https://doi.org/10.1073/pnas.88.20.9107
 Ma, S.Q., Wang, X.H., Lei, J.H. and Feng, Z.S. (2010) Dynamics of the Delay Hematological Cell Model. International Journal of Biomathematics, 3, 105-125.
 Ma, S.Q., Feng, Z.S. and Lu, Q.S. (2009) Dynamics and Double Hopf Bifurtions of the Rose-Hindmarsh Model with Time Delay. International Journal of Bifurcation and Chaos, 19, 1-19.
 Ma, S.Q., Z.S. and Lu, Q.S. (2008) A Two Parameter Criteria for Delay Differential Equations. Discrete and Continuous Dynamical Systems: Series B, 9, 397-413.
 Wang, Z., Hu, H.Y., Xu, Q. and Stepan, G. (2016) Effect of Delay Combinations on Stability and Hopf Bifurcation of an Oscillator with Acceleration-Derivative Feedback. International Journal of Nonlinear Mechanics, 94, 392-399.