Bifurcation Analysis of the Regulatory Modules of the Mammalian M/G₁ Phase with Time Delay

S. Q. Ma

Department of Mathematics, China Agricultural University, Beijing, China.

DOI: 10.4236/ijmnta.2021.102003                

Abstract

G₀/G₁ “gaps” joint the S phase and M phase to form the cell cycle. The dynamics of enzyme reaction to drive the target protein production in M phase is analyzed mathematically. Time delay is introduced since the signal transmission need time in G₀/G₁ “gaps” phase. Hopf bifurcation of DDEs model is analyzed by applying geometrical analytical method. The instability oscillating periodic solutions arise as subcritical Hopf bifurcation occurs. The Hysteresis phenomena of the limit cycle are also observed underlying the saddle-node bifurcation of the limit cycle. Due to stability switching, interestingly, the bifurcating periodical solution dies out near the vicinity of Hopf lines. By Lyapunov-Schmidt reduction scheme, the normal form is computed on the center manifold. Finally, it is verified that the theory analytical results are in coincidence with the numerical simulation.

Keywords

Cells Cycle, Time Delay, Hopf Bifurcation, Stability Switching, Normal Form

1. Introduction

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 [1] [2] [3] [4], 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 G₀/G₁ phase [5].

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 [6] [7] [8]. As is well known, G₀/G₁ “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 [9] [10]. Hence, Cdk2 dimer dynamics is dominated by the well-known Michaelis-Menten rate laws. The Michaelis-Menten rate law is described as

$\text{rate of loss of substrate}=\frac{V_{\max }S}{K_m+S}$

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 [11] [12]. Based on the above biology background, we set up delay differential equations (DDEs) to express protein regulation in the cell cycle, and write G₁ mitosis phase into a useful mathematical expression.

$\frac{\text{d}y}{\text{d}t}=-k_{1}y+k_2\frac{s_0+y{\left(t-\tau \right)}^2}{{\theta }_1+y{\left(t-\tau \right)}^2}$ (1.1)

Besides, one mother cell can be divided into two daughter cells during mitosis stage; the dynamics in G₁/M phase can also be expressed as

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=-\gamma x+g\left(x,y\right)\\ \frac{\text{d}y}{\text{d}t}=\delta -k_{1}y+2{\text{e}}^{-r\tau }g\left(x\left(t-\tau \right),y\left(t-\tau \right)\right)+k_3\left(y\left(t-{\tau }_2\right)-y\left(t\right)\right)\end{array}$ (1.2)

wherein x denotes enzyme concentration, and y represents gene protein concentration, $\tau$ is time delay with its valuable estimation between 2.8 h and 7.8 h. The associated function $g\left(x\mathrm{,}y\right)$ is written as

$\begin{array}{l}g\left(x,y\right)=k_2\frac{s_0+y^2}{{\theta }_1+y^2}\frac{{\theta }_2}{{\theta }_2+x}\hfill \end{array}$ (1.3)

State feedback control is to introduce concentration difference into the model which reflects the time delay effects on dynamical evolution behavior. Herein ${\tau }_2=\tau +{\tau }_1$, with ${\tau }_1$ denotes time delay in “gap” junction as G₀ 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 [11] [12], 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 $\tau$. 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 [13] [14] [15] [16] [17] 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 [16] [17] [18], 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 G₁ Phase of Cell Growth

The bistability phenomena is observed as varying free parameter pairs $\left(k_1\mathrm{,}{\theta }_1\right)$ 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

$\left(\Gamma \right):\{\begin{array}{l}-k_1{x}^3-k_1{\theta }_{1}x+x^2+s_0=0;\hfill \\ -3k_1{x}^2-k_1{\theta }_1+2x=0;\hfill \end{array}$ (2.1)

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 $\Gamma$ is unstable region lying on three-dimensional space, and the curve $\Gamma$ 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 $\left(k_1\mathrm{,}{\theta }_1\right)$ is attached to the vicinity boundary of fold bifurcation curves. The cusp bifurcation point CP is denoted as $\left(k_c\mathrm{,}{\theta }_c\right)$, system (1.1) exhibits bistability property if choosing free parameter $k_1$ less than $k_c$.

3. Hopf Bifurcation in G₁/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 $E^{\mathrm{<mo>\; *\; </mo>}}=\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ which satisfies

$\begin{array}{l}\left(-\gamma \left(y^{*2}+{\theta }_1\right)x^{*}+k_2\left(y^{*2}+s_0\right)\right){\theta }_2-x^{*2}\gamma \left(y^{*2}+{\theta }_1\right)=0,\\ 2{\text{e}}^{-\gamma \tau }{k}_2\left(y^{*2}+s_0\right){\theta }_2-\left(y^{*2}+{\theta }_1\right)\left(k_1{y}^{*}-\delta \right)\left({\theta }_2+x^{*}\right)=0\end{array}$ (3.1)

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 $E^{\mathrm{<mo>\; *\; </mo>}}=\left(x^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ is given as

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\left(g_x-\gamma \right)x+g_{y}y\\ \frac{\text{d}y}{\text{d}t}=2{\text{e}}^{-\gamma \tau }{g}_{x}x\left(t-\tau \right)-k_{1}y+2{\text{e}}^{-\gamma \tau }{g}_{y}y\left(t-\tau \right)+k_3\left(y\left(t-{\tau }_2\right)-y\left(t\right)\right)\end{array}$ (3.2)

Therefore, the characteristic equation of the linearized version (3.2) is written as

$\Delta \left(\lambda ,k_1,\tau \right)=\left|\begin{array}{cc}{g}_x-\gamma -\lambda & g_y\\ 2{\text{e}}^{-\left(\lambda +\gamma \right)\tau }{g}_x& -k_1-k_3-\lambda +k_3{\text{e}}^{-\lambda {\tau }_2}+2{\text{e}}^{-\left(\lambda +\gamma \right)\tau }{g}_y\end{array}\right|=0$ (3.3)

Figure 2. Hopf bifurcation of the positive equilibrium solution. (a) The continuation of equilibrium solution as varying free parameter. Hopf bifurcation occurs at $k_1=0.8504$ or $k_1=0.157$, other parameters are fixed as $k_2=0.5$, ${\theta }_2=0.7$, $\delta =0.045$, $r=0.006$, ${\theta }_1=0.6$, $s=0.029$, $\gamma =0.63$, $k_3=-0.6$, $\tau =6$, ${\tau }_2=2.309$ ; (b) The curve of angle ${\varphi }_L$ and ${\varphi }_S$ with respect to the parameter $\theta$ is plotted with $k_1=0.8504$. The intersection point is associated with Hopf bifurcation; (c) The curve of $\tau \left(\theta \mathrm{,}{\varphi }_L\right)$ and $\tau \left(\theta \mathrm{,}{\varphi }_S\right)$ respectively with respect to $\theta$ is plotted with $k_1=0.8504$ ; (d) The curve of $\tau \left(\theta \mathrm{,}{\varphi }_L\right)$ and $\tau \left(\theta \mathrm{,}{\varphi }_S\right)$ respectively with respect to $\theta$ is plotted with $k_1=0.157$.

Furthermore, one obtains that

$\Delta \left(\lambda ,k_1,\tau \right)=P\left(\lambda ,c,\tau \right){\text{e}}^{-\left(\lambda +\gamma \right)\tau }+S\left(\lambda ,c,\tau \right){\text{e}}^{-\lambda {\tau }_2}+Q\left(\lambda ,c,\tau \right)$ (3.4)

with

$\begin{array}{l}P\left(\lambda ,k_1,\tau \right)=-2\gamma g_y-2g_y\lambda ,S\left(\lambda ,c,\tau \right)=-k_3\lambda ,\\ Q\left(\lambda ,k_1,\tau \right)=\gamma k_1+\gamma \lambda -g_x{k}_1-g_x\lambda +k_1\lambda +{\lambda }^2-k_3\left(g_x-\gamma -\lambda \right).\end{array}$ (3.5)

where $k_1$ is the free parameter.

Suppose $\lambda =i\omega \left(\omega >0\right)$, substitute it into Equation (3.5) and separate the real part from the imaginary part to get

$\begin{array}{l}{P}_R{\text{e}}^{-r\tau }\cos \left(\omega \tau \right)+P_I{\text{e}}^{-r\tau }\sin \left(\omega \tau \right)+S_R\cos \left(\omega {\tau }_2\right)+S_I\sin \left(\omega {\tau }_2\right)+Q_R=0,\\ P_I{\text{e}}^{-r\tau }\cos \left(\omega \tau \right)-P_R{\text{e}}^{-r\tau }\sin \left(\omega \tau \right)-S_R\sin \left(\omega {\tau }_2\right)+S_I\cos \left(\omega {\tau }_2\right)+Q_I=0.\end{array}$ (3.6)

with $P_R\mathrm{,}{P}_I\mathrm{,}{Q}_R\mathrm{,}{Q}_I$ represents the real part and the imaginary part of the polynomial $P\mathrm{,}Q$ respectively.

Solving $\cos \left(\omega \tau \right)\mathrm{,}\sin \left(\omega \tau \right)$ from Equations (3.6) to get

$\begin{array}{l}\cos \left(\omega \tau \right)=-\frac{\left(P_I{Q}_I+P_R{Q}_R\right)+\cos \left(\omega {\tau }_2\right)\left(P_I{S}_I+P_R{S}_R\right)+\sin \left(\omega {\tau }_2\right)\left(P_R{S}_I-P_I{S}_R\right)}{{\text{e}}^{-r\tau }\left(P_I^2+P_R^2\right)},\\ \sin \left(\omega \tau \right)=-\frac{\left(P_I{Q}_R-P_R{Q}_I\right)+\cos \left(\omega {\tau }_2\right)\left(P_I{S}_R-P_R{S}_I\right)+\sin \left(\omega {\tau }_2\right)\left(P_I{S}_I+P_R{S}_R\right)}{{\text{e}}^{-r\tau }\left(P_I^2+P_R^2\right)}\end{array}$ (3.7)

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 $s$ and the associated angle variables $\theta =s\omega ,\varphi =\omega \tau ,\psi =\omega {\tau }_2$, we define the new mapping

$\theta \to {\theta }_0+2l\pi ,\varphi ={\varphi }_0+2k\pi ,\psi ={\psi }_0+2m\pi$

with $l\mathrm{,}k\mathrm{,}m\in N_0^{+}$, then rewrite the polynomial $P_R\mathrm{,}{P}_I\mathrm{,}{Q}_R\mathrm{,}{Q}_I$ as

$\begin{array}{l}{P}_R\left(\omega ,k_1,\tau \right)=P_R\left(\frac{\theta }{s},k_1,\frac{\varphi s}{\theta }\right),P_I\left(\omega ,k_1,\tau \right)=P_I\left(\frac{\theta }{s},k_1,\frac{\varphi s}{\theta }\right),\\ Q_R\left(\omega ,k_1,\tau \right)=Q_R\left(\frac{\theta }{s},k_1,\frac{\varphi s}{\theta }\right),Q_I\left(\omega ,k_1,\tau \right)=Q_I\left(\frac{\theta }{s},k_1,\frac{\varphi s}{\theta }\right),\\ S_R\left(\omega ,k_1,\tau \right)=S_R\left(\frac{\theta }{s},k_1,\frac{\varphi s}{\theta }\right),S_I\left(\omega ,c,\tau \right)=S_I\left(\frac{\theta }{s},k_1,\frac{\varphi s}{\theta }\right)\end{array}$ (3.9)

Define the new function

$\begin{array}{l}F\left(\theta ,k_1,s\right)={\left(P_I{Q}_I+P_R{Q}_R\right)}^2+{\left(P_I{Q}_R-P_R{Q}_I\right)}^2-{\text{e}}^{-2r\frac{\varphi s}{\theta }}{\left(P_I^2+P_R^2\right)}^2\\ +2\left(P_I{Q}_I+P_R{Q}_R\right)\left(P_I{S}_I+P_R{S}_R\right)\cos \left(\frac{\theta {\tau }_2}{s}\right)\\ +2\left(P_I{Q}_R-P_R{Q}_I\right)\left(P_I{S}_R-P_R{S}_I\right)\cos \left(\frac{\theta {\tau }_2}{s}\right)\end{array}$ (3.10)

$\begin{array}{l}+2\left(P_I{Q}_I+P_R{Q}_R\right)\left(P_R{S}_I-P_I{S}_R\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)\\ +2\left(P_I{Q}_R-P_R{Q}_I\right)\left(P_I{S}_I+P_R{S}_R\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)\\ +{\left(P_I{S}_I+P_R{S}_R\right)}^2+{\left(P_I{S}_R-P_R{S}_I\right)}^2\end{array}$

By Equation (3.10), one defines curve L by

$L:\varphi =-\frac{\theta }{2rs}\ln \frac{H\left(P_R,P_I,Q_R,Q_I,S_R,S_I,\cos \left(\frac{\theta {\tau }_2}{s}\right),\sin \left(\frac{\theta {\tau }_2}{s}\right)\right)}{{\left(P_I^2+P_R^2\right)}^2}$ (3.11)

with

$\begin{array}{l}H\left(P_R,P_I,Q_R,Q_I,S_R,S_I,\cos \left(\frac{\theta {\tau }_2}{s}\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)\right)\\ ={\left(P_I{Q}_I+P_R{Q}_R\right)}^2+{\left(P_I{Q}_R-P_R{Q}_I\right)}^2+2\left(P_I{Q}_I+P_R{Q}_R\right)\left(P_I{S}_I+P_R{S}_R\right)\cos \left(\frac{\theta {\tau }_2}{s}\right)\\ +2\left(P_I{Q}_R-P_R{Q}_I\right)\left(P_I{S}_R-P_R{S}_I\right)\cos \left(\frac{\theta {\tau }_2}{s}\right)\\ +2\left(P_I{Q}_I+P_R{Q}_R\right)\left(P_R{S}_I-P_I{S}_R\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)\\ +2\left(P_I{Q}_R-P_R{Q}_I\right)\left(P_I{S}_I+P_R{S}_R\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)+{\left(P_I{S}_I+P_R{S}_R\right)}^2+{\left(P_I{S}_R-P_R{S}_I\right)}^2\end{array}$

By Equation (3.11), we also get the description of curve S

$S:\varphi =\arctan \frac{U\left(P_R,P_I,Q_R,Q_I,S_R,S_I,\cos \left(\frac{\theta {\tau }_2}{s}\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)\right)}{V\left(P_R,P_I,Q_R,Q_I,S_R,S_I,\cos \left(\frac{\theta {\tau }_2}{s}\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)\right)}-2k\pi$ (3.12)

with

$\begin{array}{l}U\left(P_R,P_I,Q_R,Q_I,S_R,S_I,\cos \left(\frac{\theta {\tau }_2}{s}\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)\right)\\ =\left(P_I{Q}_R-P_R{Q}_I\right)+\cos \left(\frac{\theta {\tau }_2}{s}\right)\left(P_I{S}_R-P_R{S}_I\right)+\sin \left(\frac{\theta {\tau }_2}{s}\right)\left(P_I{S}_I+P_R{S}_R\right),\\ V\left(P_R,P_I,Q_R,Q_I,S_R,S_I,\cos \left(\frac{\theta {\tau }_2}{s}\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)\right)\\ =\left(P_I{Q}_I+P_R{Q}_R\right)+\cos \left(\frac{\theta {\tau }_2}{s}\right)\left(P_I{S}_I+P_R{S}_R\right)+\sin \left(\frac{\theta {\tau }_2}{s}\right)\left(P_R{S}_I-P_I{S}_R\right).\end{array}$

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 ${\varphi }_L={\varphi }_0+2l\pi$ and ${\varphi }_S+k\pi ={\varphi }_L$ for some k and l. With the assumption that the intersection of two curves is expressed as $\left({\theta }^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\varphi }^{\mathrm{<mo>\; *\; </mo>}}\right)$ which satisfies ${\varphi }_S\left({\theta }^{\mathrm{<mo>\; *\; </mo>}}\right)+k\pi ={\varphi }_L\left({\theta }^{\mathrm{<mo>\; *\; </mo>}}\right)$, we compute the time delay ${\tau }^{\mathrm{<mo>\; *\; </mo>}}$ as the critical value of Hopf bifurcation, ${\tau }^{*}=\frac{{\varphi }_{L}s}{\theta }$. Note that $\tau$ is represented by the introduction delay s. To calculate the Hopf bifurcation line, Equation (3.11) builts the inverse function $c=c\left(\theta \right)$ with ${\theta }_0$ lying inside some subinterval of $\left[\mathrm{0,2}\pi \right]$ since we have defined the relationship between $\theta$ and ${\theta }_0$ by mapping (3.8), and Hopf bifurcation occurs at

$\tau \left(\theta \right)=\frac{s}{\theta }\left(\arctan \frac{U\left(P_R,P_I,Q_R,Q_I,S_R,S_I,\cos \left(\frac{\theta {\tau }_2}{s}\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)\right)}{V\left(P_R,P_I,Q_R,Q_I,S_R,S_I,\cos \left(\frac{\theta {\tau }_2}{s}\right)\sin \left(\frac{\theta {\tau }_2}{s}\right)\right)}-2k\pi \right)$ (3.13)

Following we compute the transversal condition of Hopf bifurcation, differentiate the rightside of Equation (3.4) with respect to $\tau$, one obtains

$\begin{array}{l}\left(\frac{\partial P}{\partial \lambda }\frac{\text{d}\lambda }{\text{d}\tau }+\frac{\partial P}{\partial \tau }\right){\text{e}}^{-\left(\lambda +r\right)\tau }+P{\text{e}}^{-\left(\lambda +r\right)\tau }\left(-\frac{\text{d}\lambda }{\text{d}\tau }\tau -r-\lambda \right)\\ +\left(\frac{\partial S}{\partial \lambda }\frac{\text{d}\lambda }{\text{d}\tau }+\frac{\partial S}{\partial \tau }-{\tau }_{2}S\frac{\text{d}\lambda }{\text{d}\tau }\right){\text{e}}^{-\lambda {\tau }_2}+\left(\frac{\partial Q}{\partial \lambda }\frac{\text{d}\lambda }{\text{d}\tau }+\frac{\partial Q}{\partial \tau }\right)=0\end{array}$ (3.14)

then by solving $\frac{\text{d}\lambda }{\text{d}\tau }$ from Equation (3.14), we get

$\frac{1}{\frac{\text{d}\lambda }{\text{d}\tau }}=\frac{\frac{\partial P}{\partial \lambda }{\text{e}}^{-\left(\lambda +r\right)\tau }-P{\text{e}}^{-\left(\lambda +r\right)\tau }+\left(\frac{\partial S}{\partial \lambda }-{\tau }_{2}S\right){\text{e}}^{-\lambda {\tau }_2}+\frac{\partial Q}{\partial \lambda }}{\frac{\partial P}{\partial \tau }{\text{e}}^{-\left(\lambda +r\right)\tau }+P{\text{e}}^{-\left(\lambda +r\right)\tau }\left(-r-\lambda \right)+\frac{\partial S}{\partial \tau }{\text{e}}^{-\lambda {\tau }_2}+\frac{\partial Q}{\partial \tau }}$ (3.15)

differentiate the rightside of Equation (3.4) with respect to $\theta$ to further get

$\begin{array}{l}\left(\frac{\partial P}{\partial \lambda }\frac{\text{d}\lambda }{\text{d}\theta }+\frac{\partial P}{\partial c}\frac{\text{d}c}{\text{d}\theta }+\frac{\partial P}{\partial \tau }\frac{\text{d}\tau }{\text{d}\theta }\right){\text{e}}^{-\left(\lambda +r\right)\tau }+P{\text{e}}^{-\left(\lambda +r\right)\tau }\left(-\frac{\text{d}\lambda }{\text{d}\theta }\tau -\left(r+\lambda \right)\frac{\text{d}\tau }{\text{d}\theta }\right)\\ -{\tau }_{2}S\frac{\text{d}\lambda }{\text{d}\theta }{\text{e}}^{-\lambda {\tau }_2}+\left(\frac{\partial S}{\partial \lambda }\frac{\text{d}\lambda }{\text{d}\theta }+\frac{\partial S}{\partial c}\frac{\text{d}c}{\text{d}\theta }+\frac{\partial S}{\partial \tau }\frac{\text{d}\tau }{\text{d}\theta }\right){\text{e}}^{-\lambda {\tau }_2}\\ +\left(\frac{\partial Q}{\partial \lambda }\frac{\text{d}\lambda }{\text{d}\theta }+\frac{\partial Q}{\partial c}\frac{\text{d}c}{\text{d}\theta }+\frac{\partial Q}{\partial \tau }\frac{\text{d}\tau }{\text{d}\theta }\right)=0\end{array}$ (3.16)

Noticed that $\lambda =i\omega$ and $\omega s=\theta$, we get $\frac{\text{d}\lambda }{\text{d}\theta }=i\frac{1}{s}$, then substitute it into Equation (3.16) to get

$\begin{array}{l}\left(\frac{\partial P}{\partial \lambda }\frac{1}{s}i+\frac{\partial P}{\partial c}\frac{\text{d}c}{\text{d}\theta }\right){\text{e}}^{-\left(\lambda +r\right)\tau }+P{\text{e}}^{-\left(\lambda +r\right)\tau }\left(-i\frac{1}{s}\tau \right)+S{\text{e}}^{-\lambda {\tau }_2}\left(-i\frac{1}{s}{\tau }_2\right)\\ +\left(\frac{\partial S}{\partial \lambda }\frac{1}{s}i+\frac{\partial S}{\partial c}\frac{\text{d}c}{\text{d}\theta }\right){\text{e}}^{-\lambda {\tau }_2}+\left(\frac{\partial Q}{\partial \lambda }\frac{1}{s}i+\frac{\partial Q}{\partial c}\frac{\text{d}c}{\text{d}\theta }\right)\\ =\left(-\frac{\partial P}{\partial \tau }{\text{e}}^{-\left(\lambda +r\right)\tau }+P{\text{e}}^{-\left(\lambda +r\right)\tau }\left(r+\lambda \right)-\frac{\partial S}{\partial \tau }{\text{e}}^{-\lambda {\tau }_2}-\frac{\partial Q}{\partial \tau }\right)\frac{\text{d}\tau }{\text{d}\theta }\end{array}$ (3.17)

then by Equation (3.17) and Equation (3.15), one has

$\begin{array}{l}\Re \left({\frac{1}{\frac{\text{d}\lambda }{\text{d}\tau }}|}_{\lambda =i\omega }\right)=-\Re (\left(\frac{\partial P}{\partial \lambda }{\text{e}}^{-\left(\lambda +r\right)\tau }-P\tau {\text{e}}^{-\left(\lambda +r\right)\tau }-{\tau }_{2}S{\text{e}}^{-\lambda {\tau }_2}+\frac{\partial S}{\partial \lambda }{\text{e}}^{-\lambda {\tau }_2}+\frac{\partial Q}{\partial \lambda }\right)\\ \left(abcc.\frac{1}{s}\left(-i\right)+sbcc.\frac{\text{d}c}{\text{d}\theta }\right))\end{array}$ (3.18)

herein $abcc\mathrm{.}$ denotes the conjugate part of $\frac{\partial P}{\partial \lambda }{\text{e}}^{-\left(\lambda +r\right)\tau }-P\tau {\text{e}}^{-\left(\lambda +r\right)\tau }-{\tau }_{2}S{\text{e}}^{-\lambda {\tau }_2}+\frac{\partial S}{\partial \lambda }{\text{e}}^{-\lambda {\tau }_2}+\frac{\partial Q}{\partial \lambda }$, and $sbcc\mathrm{.}$ expresses the conjugate part of $\frac{\partial P}{\partial c}{\text{e}}^{-\left(\lambda +r\right)\tau }+\frac{\partial S}{\partial c}{\text{e}}^{-\lambda {\tau }_2}+\frac{\partial Q}{\partial c}$. Noticed the expression of the characteristic polynomial (3.4), one has

$\Re \left({\frac{1}{\frac{\text{d}\lambda }{\text{d}\tau }}|}_{\lambda =i\omega }\right)=-sign\Re \left({\Delta }_{\lambda }\left(\lambda ,c,\tau \right){\bar{\Delta }}_c\left(\lambda ,c,\tau \right)\right)sign\left(\frac{\text{d}c}{\text{d}\theta }\right)$ (3.19)

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 $k_1$ 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 $\left(k_3\mathrm{,}\tau \right)$ plane. (a) The stability switching phenomena for the equilibrium solution with chosen $k_1=0.1577$ ; (b) Hopf line on parameter $\left(k_3\mathrm{,}\tau \right)$ plane with $k_1=0.1557$ ; (c) The stability switching for the equilibrium solution as $k_1=0.85$ ; (d) Hopf line on parameter $\left(k_3\mathrm{,}\tau \right)$ plane with $k_1=0.85$.

S satisfies ${\Phi }_L=\mathrm{mod}\left({\Phi }_S\mathrm{,2}\pi \right)$, as shown in Figure 2(b). Hopf bifurcation occurs which specifies angle parameter $\theta$ to derive time delay $\tau$. The curve S and L are drawn in $\left(\theta \mathrm{,}\varphi \right)$ plane. By equality $\tau \triangleq \tau \left(\theta ,{\varphi }_L\right)=\tau \triangleq \tau \left(\theta ,{\varphi }_S\right)$, 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 $\tau$, as shown in Figure 3. In Figure 3(b) and Figure 3(d), Hopf lines are drawn on $\left(k_3-\tau \right)$ 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 $k_1$, Hopf bifurcation occurs at $k_1=0.876$ and $k_1=0.157$, 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 $k_1=0.876$ ; (b) near $k_1=0.157$ ; (c) near $\tau =21.8$ ; (d) near $\tau =31.76$.

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 $k_1$ 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 $\left(k_3,\tau \right)=\left(k_3^{*},{\tau }_c\right)$ 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 $0<\epsilon \ll 1$, set $k_3-k_{3c}=\epsilon k_{\epsilon }$, $\tau -{\tau }_c=\epsilon {\tau }_{\epsilon }$ and scale $x-x^{*}=\sqrt{\epsilon }x$, $y-y^{*}=\sqrt{\epsilon }y$, then system (1.2) is transformed into its 3rd Taylor expansion truncation

$\begin{array}{l}\frac{\text{d}x}{\text{d}t}=-\gamma x+g_{x}x+g_{y}y+\sqrt{\epsilon }\frac{1}{2}{g}_{xx}{x}^2+\sqrt{\epsilon }{g}_{xy}xy+\sqrt{\epsilon }\frac{1}{2}{g}_{yy}{y}^2\\ +\epsilon \frac{1}{6}{g}_{xxx}{x}^3+\epsilon \frac{1}{2}{g}_{xxy}{x}^{2}y+\epsilon \frac{1}{2}{g}_{xyy}x{y}^2+\epsilon \frac{1}{6}{g}_{yyy}{y}^3\end{array}$

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 $k_1=0.1574,\tau =16$ ; (b) The corresponding phase portraits of the period-2 solution; (c) The time series solution of the quasi-periodical solutions with $k_1=0.1574,\tau =23$ ; (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 $k_1=0.85,\tau =16$ ; (b) The corresponding phase portraits of the period-2 solution; (c) The time series solution of the quasi-periodical solutions with $k_1=0.85,\tau =23$ ; (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 $k_1=0.1574$ ; (b) The Poincare section with $k_1=0.85$.

$\frac{\text{d}y}{\text{d}t}=2{\text{e}}^{-r{\tau }_c}{g}_{x}x\left(t-{\tau }_c\right)-2r\epsilon {\tau }_{\epsilon }{\text{e}}^{-r{\tau }_c}{g}_{x}x\left(t-{\tau }_c\right)-\left(k_1+k_{3c}\right)y$

$\begin{array}{l}+2{\text{e}}^{-r{\tau }_c}{g}_x\left(x\left(t-{\tau }_c-\epsilon {\tau }_{\epsilon }\right)-x\left(t-{\tau }_c\right)\right)+k_{3}y\left(t-{\tau }_2\right)\\ +\epsilon k_{\epsilon }y\left(t-{\tau }_2\right)+2{\text{e}}^{-r{\tau }_c}{g}_{y}y\left(t-{\tau }_c\right)-2r\epsilon {\tau }_{\epsilon }{g}_{y}y\left(t-{\tau }_c\right)\\ +2{\text{e}}^{-r{\tau }_c}{g}_y\left(y\left(t-{\tau }_c-\epsilon {\tau }_{\epsilon }\right)-y\left(t-{\tau }_c\right)\right)\end{array}$

$\begin{array}{l}+\sqrt{\epsilon }{\text{e}}^{-r{\tau }_c}{g}_{xx}x{\left(t-{\tau }_c\right)}^2+\sqrt{\epsilon }{\text{e}}^{-r{\tau }_c}{g}_{xy}x\left(t-{\tau }_c\right)y\left(t-{\tau }_c\right)\\ +\sqrt{\epsilon }{\text{e}}^{-r{\tau }_c}{g}_{yy}y{\left(t-{\tau }_c\right)}^2+\epsilon \frac{1}{3}{\text{e}}^{-r{\tau }_c}{g}_{xxx}x{\left(t-{\tau }_c\right)}^3\\ +\epsilon {\text{e}}^{-r{\tau }_c}{g}_{xxy}x{\left(t-{\tau }_c\right)}^{2}y\left(t-{\tau }_c\right)+\epsilon {\text{e}}^{-r{\tau }_c}{g}_{xyy}x\left(t-{\tau }_c\right)y{\left(t-{\tau }_c\right)}^2\\ +\epsilon \frac{1}{3}{\text{e}}^{-r{\tau }_c}{g}_{yyy}y{\left(t-{\tau }_c\right)}^3\end{array}$ (4.1)

with initial value definition

$\begin{array}{l}\left(\begin{array}{c}x\left(t\right)\\ y\left(t\right)\end{array}\right)=\left(\begin{array}{c}{\varphi }_1\\ {\varphi }_2\end{array}\right),-{\tau }_m\le t\le 0\hfill \end{array}$

wherein ${\tau }_m=\max \left({\tau }_2,\tau \right)$. According to the property of DDEs, solutions of Equation (4.1) are continuous on Banach space $C=\left(\left[-{\tau }_m,0\right],R^2\right)$. With few discontinuity jumps, the solution operator is differentiable on the extended phase space $BC=\left(\left[-{\tau }_m,0\right],R^2\right)$. With the assumption $Z={\left(x\mathrm{,}y\right)}^{\text{T}}$, and denote $Z\left(t+s\right)=Z_t\left(s\right)$ for any $s\in \left[-{\tau }_m\mathrm{,0}\right]$, then Equation (4.1) is written as

$\frac{\text{d}Z\left(t\right)}{\text{d}t}=L\left(\epsilon \right)Z_t+F\left(Z_t\left(0\right)\mathrm{,}{Z}_t\left(-\tau \right)\mathrm{,}{z}_t\left(-{\tau }_2\right)\mathrm{,}\epsilon \right)$ (4.2)

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 $\varphi \in BC$,

$L\left(\epsilon \right)\varphi ={\displaystyle {\int }_{-{\tau }_c}^0}\text{d}\eta \left(\theta \right)\varphi +\epsilon {\displaystyle {\int }_{-{\tau }_m}^0}\text{d}{\eta }_1\left(\theta \right)\varphi +{\displaystyle {\int }_{-{\tau }_m}^0}\text{d}{\eta }_2\left(\theta \right)\varphi ,$ (4.3)

with

$\text{d}\eta \left(\theta \right)=\left(\begin{array}{cc}\left(-\gamma +g_x\right)\delta \left(\theta \right)& g_y\delta \left(\theta \right)\\ 2{\text{e}}^{-r{\tau }_c}{g}_x\delta \left(\theta +{\tau }_c\right)& -\left(k_1+k_{3c}\right)\delta \left(\theta \right)+k_3\delta \left(\theta +{\tau }_2\right)+2{\text{e}}^{-r{\tau }_c}{g}_y\delta \left(\theta +{\tau }_c\right)\end{array}\right)$

$\text{d}{\eta }_1\left(\theta \right)=\left(\begin{array}{cc}0& 0\\ -2r{\tau }_{\epsilon }{\text{e}}^{-r{\tau }_c}{g}_x\delta \left(\theta +{\tau }_c\right)& k_{\epsilon }\delta \left(\theta +{\tau }_2\right)-2r{\tau }_{\epsilon }{g}_y\delta \left(\theta +{\tau }_c\right)\end{array}\right)$

$\text{d}{\eta }_2\left(\theta \right)=\left(\begin{array}{cc}0& 0\\ 2{\text{e}}^{-r{\tau }_c}{g}_x\left(\delta \left(\theta +{\tau }_c+\epsilon {\tau }_{\epsilon }\right)-\delta \left(\theta +{\tau }_c\right)\right)& 2{\text{e}}^{-r{\tau }_c}{g}_y\left(\delta \left(\theta +{\tau }_c+\epsilon {\tau }_{\epsilon }\right)-\delta \left(\theta +{\tau }_c\right)\right)\end{array}\right)$

and the nonlinear part $F(\cdot )$ is written as

$\begin{array}{c}F\left(\varphi \left(0\right),\varphi \left(-\tau \right)\right)=(\frac{1}{2}{g}_{xx}{\varphi }_1{\left(0\right)}^2+g_{xy}{\varphi }_1\left(0\right){\varphi }_2\left(0\right)+\frac{1}{2}{g}_{yy}{\varphi }_2{\left(0\right)}^2+\frac{1}{6}{g}_{xxx}{\varphi }_1{\left(0\right)}^3\\ +\frac{1}{2}{g}_{xxy}{\varphi }_1{\left(0\right)}^2{\varphi }_2\left(0\right)+\frac{1}{2}{g}_{xyy}{\varphi }_1\left(0\right){\varphi }_2{\left(0\right)}^2+\frac{1}{6}{g}_{yyy}{\varphi }_2{\left(0\right)}^3,\\ {\text{e}}^{-r{\tau }_c}{g}_{xx}{\varphi }_1{\left(-{\tau }_c\right)}^2+{\text{e}}^{-r{\tau }_c}{g}_{xy}{\varphi }_1\left(-{\tau }_c\right){\varphi }_2\left(-{\tau }_c\right)+{\text{e}}^{-r{\tau }_c}{g}_{yy}{\varphi }_2{\left(-{\tau }_c\right)}^2\end{array}$ $\begin{array}{c}+\frac{1}{3}{\text{e}}^{-r{\tau }_c}{g}_{xxx}{\varphi }_1{\left(-{\tau }_c\right)}^3+{\text{e}}^{-r{\tau }_c}{g}_{xxy}{\varphi }_1{\left(-{\tau }_c\right)}^2{\varphi }_2\left(-{\tau }_c\right)\\ +{\text{e}}^{-r{\tau }_c}{g}_{xyy}{\varphi }_1\left(-{\tau }_c\right){\varphi }_2{\left(-{\tau }_c\right)}^2+\frac{1}{3}{\text{e}}^{-r{\tau }_c}{g}_{yyy}{\varphi }_2{\left(-{\tau }_c\right)}^3)\end{array}$ (4.4)

Consider the linear operator (4.3), it’s an infinitesimal generator of the strong continuous semigroup in phase space BC, and for any $\varphi \in BC$, we define new linear operator $A\left(\epsilon \right)\mathrm{<mo>\; :\; </mo>}BC\to BC$ and its adjoint operator $A^{\mathrm{<mo>\; *\; </mo>}}\left(\epsilon \right)\mathrm{<mo>\; :\; </mo>}B{C}^{\mathrm{<mo>\; *\; </mo>}}\to B{C}^{\mathrm{<mo>\; *\; </mo>}}$ wherein $B{C}^{\mathrm{<mo>\; *\; </mo>}}=C^{\mathrm{<mo>\; *\; </mo>}}\left(\left[\mathrm{0,}\tau \right]\mathrm{,}{R}^2\right)$, that is,

$A\left(\epsilon \right)\varphi =\{\begin{array}{l}\frac{\text{d}\varphi }{\text{d}\theta },\text{for}-\tau \le \theta <0\hfill \\ L\left(\epsilon \right)\varphi ,\text{for}\theta =0\hfill \end{array}$ (4.5)

and

$A^{*}\left(\epsilon \right)\psi =\{\begin{array}{l}-\frac{\text{d}\psi }{\text{d}s},\text{for}0<s\le \tau \\ L^{*}\left(\epsilon \right)\psi ,\text{for}s=0\end{array}$ (4.6)

with

$L^{*}\left(\epsilon \right)\psi ={\displaystyle {\int }_{-\tau }^0}\text{d}{\eta }^{\text{T}}\left(s\right)\psi \left(-s\right)$ (4.7)

For $\varphi \in BC\mathrm{,}\psi \in B{C}^{\mathrm{<mo>\; *\; </mo>}}$, define the bilinear form

$\langle \psi ,\varphi \rangle ={\bar{\psi }}^{\text{T}}\left(0\right)\varphi \left(0\right)-{\displaystyle {\int }_{-{\tau }_m}^0}{\displaystyle {\int }_0^{\theta }}{\bar{\psi }}^{\text{T}}\left(\xi -\theta \right)\text{d}\eta \left(\theta \right)\varphi \left(\xi \right)\text{d}\xi$ (4.8)

Based on the above analysis, Equation (4.2) can be written as its operator differential equations, that is

${Z^{\prime }}_t=A\left(\epsilon \right)Z_t+R{Z}_t$ (4.8)

with

$R\varphi =\{\begin{array}{l}0,\text{for}-\tau \le \theta <0\hfill \\ F\left(\varphi \left(0\right),\varphi \left(-\tau \right),\epsilon \right),\text{for}\theta =0\hfill \end{array}$ (4.9)

Considering the linear version of differential operator (4.8), Hopf bifurcation occurs at $\epsilon =0$, and the associated characteristic roots set is denoted as $\Lambda =\left\{i\omega \mathrm{,}-i\omega \right\}$ 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 $\Lambda$ and its complementary subspace Q. We suppose eigenspace P is spanned by $P=span\left\{q\left(\theta \right),\bar{q}\left(\theta \right)\right\},-\tau \le \theta \le 0$ where $q\left(\theta \right)$ and its conjugation vector $\bar{q}\left(\theta \right)$ are respectively the eigenvector of $i\omega$ and $-i\omega$. We represent the eigenbasis $\Phi \left(\theta \right)\mathrm{,}\Psi \left(s\right)$ of $\Lambda$ associated with linear operator $A\left(0\right)$ and its adjoint operator $A^{\mathrm{<mo>\; *\; </mo>}}\left(0\right)$ respectively, then

$A\left(0\right)\Phi \left(\theta \right)=\left(\begin{array}{cc}i\omega & 0\\ 0& -i\omega \end{array}\right)\Phi \left(\theta \right),A^{*}\left(0\right){\Psi }^{\text{T}}\left(s\right)=\left(\begin{array}{cc}-i\omega & 0\\ 0& i\omega \end{array}\right){\Psi }^{\text{T}}\left(s\right)$ (4.10)

By the definition of Equation (4.8), we have $\langle \Psi \mathrm{,}\Phi \rangle =I$.

With a possible discontinuity jump at $\theta =0$, we also define the mapping on the phase space as the following

$X_0=\{\begin{array}{l}I\text{for}\theta =0,\hfill \\ 0\text{for}-\tau \le \theta <0\hfill \end{array}$ (4.11)

for any $Z_t\in BC$, we can write it as $Z_t=\varphi +X_0\alpha$, define the projection mapping $\Pi \mathrm{<mo>\; :\; </mo>}BC\to P$, then we have $\Pi Z_t=\Phi \left(\theta \right)\left[\langle \Psi \mathrm{,}\varphi \rangle +{\bar{\Psi }}^{\text{T}}\left(0\right)w\right]$, therefore, we can write $Z_t$ as the expression $Z_t=\Phi \left(\begin{array}{c}z\\ \bar{z}\end{array}\right)+Y_t$ with $Y_t\in Q$. Further computation express the differentiate relationship of axis variable $z\left(t\right)$ given as

$\begin{array}{l}\left(\begin{array}{c}{z}^{\prime }\left(t\right)\\ {\bar{z}}^{\prime }\left(t\right)\end{array}\right)=\left(\begin{array}{cc}i\omega & 0\\ 0& -i\omega \end{array}\right)\left(\begin{array}{c}z\\ \bar{z}\end{array}\right)+{\bar{\Psi }}^{\text{T}}\left(0\right)\left(A\left(\epsilon \right)-A\left(0\right)\right)\Phi \left(\begin{array}{c}z\\ \bar{z}\end{array}\right)\\ +{\bar{\Psi }}^{\text{T}}\left(0\right)R\left(\Phi \left(\begin{array}{c}z\\ \bar{z}\end{array}\right)+Y_t\right),\end{array}$

$\begin{array}{l}{Y^{\prime }}_t=A\left(0\right)Y_t+\left(I-\Pi \right)X_{0}R\left(\Phi \left(\begin{array}{c}z\\ \bar{z}\end{array}\right)+Y_t\right)\\ =A\left(0\right)Y_t+\{\begin{array}{l}-\Phi \left(\theta \right){\bar{\Psi }}^{\text{T}}\left(0\right)F\left(\Phi \left(0\right)z+Y_t\left(0\right),\Phi \left(-\tau \right)z+Y_t\left(-\tau \right)\right),\text{for}-\tau \le \theta <0\hfill \\ F\left(\Phi \left(0\right)z+Y_t\left(0\right),\Phi \left(-\tau \right)z+Y_t\left(-\tau \right)\right)\hfill \\ -\Phi \left(0\right){\bar{\Psi }}^{\text{T}}\left(0\right)F\left(\Phi \left(0\right)z+Y_t\left(0\right),\Phi \left(-\tau \right)z+Y_t\left(-\tau \right)\right),\text{for}\theta =0\hfill \end{array}\end{array}$ (4.12)

Set $B=\left(\begin{array}{cc}i\omega & 0\\ 0& -i\omega \end{array}\right)$, then Equation (4.12) is written into its Taylor expansion with 3rd truncation to be expressed as

$\begin{array}{l}\left(\begin{array}{c}{z}^{\prime }\left(t\right)\\ {\bar{z}}^{\prime }\left(t\right)\end{array}\right)=B\left(\begin{array}{c}z\\ \bar{z}\end{array}\right)+f_2^{\left(11\right)}\left(z,\bar{z},0,0,\epsilon \right)+f_2^{\left(12\right)}\left(z,\bar{z},0,0,0\right)\\ +f_3^{\left(1\right)}\left(z,\bar{z},y\left(0\right),y\left(-{\tau }_c\right),0\right),\\ y^{\prime }\left(t\right)=A\left(0\right)y+\{\begin{array}{l}-\Phi \left(\theta \right)f_2^{\left(12\right)}\left(z,\bar{z},0,0,0\right),\text{for}-{\tau }_c\le \theta <0\hfill \\ f_2^{\left(2\right)}\left(z,\bar{z},0,0,0\right)-\Phi \left(0\right)f_2^{\left(12\right)}\left(z,\bar{z},0,0,0\right)\text{for}\theta =0\hfill \end{array}\end{array}$ (4.13)

and

$f_2^{\cdot }\left(z,\bar{z},0,0,\epsilon \right)={\displaystyle \underset{i+j+k=2}{\sum }}{f}_{ij0k}^{\cdot }{z}^i{\bar{z}}^j{\epsilon }^k,f_3^{\cdot }\left(z,\bar{z},0,0,0\right)={\displaystyle \underset{i+j=3}{\sum }}{f}_{ij00}^{\cdot }{z}^i{\bar{z}}^j$

Define the operator $M_2^{\mathrm{11,12,2}}$ on the homogeneous space $H_2\left(z\mathrm{,}\bar{z}\right)$ as

$\begin{array}{l}{M}_2^{1\cdot }\left(\begin{array}{c}{p}^{1\cdot }\left(z,\bar{z}\right)\\ {\bar{p}}^{1\cdot }\left(z,\bar{z}\right)\end{array}\right)=B\left(\begin{array}{c}{p}^{1\cdot }\left(z,\bar{z}\right)\\ {\bar{p}}^{1\cdot }\left(z,\bar{z}\right)\end{array}\right)-\left(\begin{array}{cc}{p}_z^{1\cdot }& p_{\bar{z}}^{1\cdot }\\ {\bar{p}}_z^{1\cdot }& {\bar{p}}_{\bar{z}}^{1\cdot }\end{array}\right)B\left(\begin{array}{c}z\\ \bar{z}\end{array}\right),\\ M_2^{2}U\left(z,\bar{z}\right)=A\left(0\right)U\left(z,\bar{z}\right)-\left(U_z,U_{\bar{z}}\right)B\left(\begin{array}{c}z\\ \bar{z}\end{array}\right)\end{array}$ (4.14)

Then the normal form on the center manifold can be written as

$\left(\begin{array}{c}{z}^{\prime }\left(t\right)\\ {\bar{z}}^{\prime }\left(t\right)\end{array}\right)=Bz+g_2^{\left(1\right)}\left(z,\bar{z},0,0,\epsilon \right)+g_3^{\left(1\right)}\left(z,\bar{z},0,0,0\right)$ (4.15)

wherein

$\begin{array}{l}{g}_2^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0,}\epsilon \right)=Pro{j}_{Im{M}_2^{11}{}^c}{f}_2^{11}\left(z\mathrm{,}\bar{z}\mathrm{,0,0,}\epsilon \right)\\ g_3^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0,0}\right)=Pro{j}_{Im{M}_3^1{}^c}{\tilde{f}}_3^{\left(1\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0,0}\right)\end{array}$ (4.16)

with

$\begin{array}{c}{\tilde{f}}_3^{\left(1\right)}\left(z,\bar{z},0,0,0\right)=f_3^{\left(1\right)}\left(z,\bar{z},0,0,0\right)-\left(\begin{array}{cc}{p}_z^{12}& p_{\bar{z}}^{12}\\ {\bar{p}}_z^{12}& {\bar{p}}_{\bar{z}}^{12}\end{array}\right)f_2^{\left(1\right)}\left(z,\bar{z},0,0,0\right)\\ +f_{2y_0}^{\left(12\right)}\left(z,\bar{z},0,0,0\right)U\left(z,\bar{z}\right)\left(0\right)\\ +f_{2y_{\tau }}^{\left(12\right)}\left(z,\bar{z},0,0,0\right)U\left(z,\bar{z}\right)\left(-{\tau }_c\right)\end{array}$ (4.17)

Hence by (4.15), one obtains the formal norm near Hopf point

$z^{\prime }\left(t\right)=i\omega z+\epsilon k_{10}z+k_{12}{z}^2\bar{z}$ (4.18)

The bases of projection $Im{M}_2^{11}\mathrm{,}Im{M}_2^{12}$ are easily to compute as

$\left(\begin{array}{c}2i\omega \bar{z}\epsilon \\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -2i\omega z\epsilon \end{array}\right)\mathrm{,}\left(\begin{array}{c}-i\omega z^2\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}i\omega z\bar{z}\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}3i\omega {\bar{z}}^2\\ 0\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -3i\omega z^2\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ -i\omega z\bar{z}\end{array}\right)\mathrm{,}\left(\begin{array}{c}0\\ i\omega {\bar{z}}^2\end{array}\right)$

The bases of $Im\left(M_1^3\right)$ are expressed as

$\begin{array}{l}\left(\begin{array}{c}-2i\omega z^3\\ 0\end{array}\right),\left(\begin{array}{c}2i\omega z{\bar{z}}^2\\ 0\end{array}\right),\left(\begin{array}{c}4i\omega {\bar{z}}^3\\ 0\end{array}\right),\left(\begin{array}{c}0\\ -4i\omega z^3\end{array}\right),\left(\begin{array}{c}0\\ -2i\omega z^2\bar{z}\end{array}\right),\left(\begin{array}{c}0\\ 2i\omega {\bar{z}}^3\end{array}\right)\hfill \end{array}$

Hence, we get

$k_{10}=f_{10001,1}^{\left(11\right)}\left(z,\bar{z},0,0\right)$ (4.19)

By choosing

$\left(\begin{array}{c}{p}^{\left(12\right)}\left(z,\bar{z}\right)\\ {\bar{p}}^{\left(12\right)}\left(z,\bar{z}\right)\end{array}\right)=\frac{f_{20000}^{\left(12\right)}}{-i\omega }{z}^2+\frac{f_{11000}^{\left(12\right)}}{i\omega }z\bar{z}+\frac{f_{02000}^{\left(12\right)}}{3i\omega }{\bar{z}}^2,$ (4.20)

we eliminate the term $f_2^{\left(12\right)}\left(z\mathrm{,}\bar{z}\mathrm{,0,0,0}\right)$. Set

$U\left(z,\bar{z}\right)=H_{2000}{z}^2+H_{1100}z\bar{z}+H_{0200}{\bar{z}}^2$ (4.21)

then by Equation (4.13) and Equation (4.14), for $-\tau \le \theta <0$, we have

$\begin{array}{l}A\left(0\right)H_{2000}\left(\theta \right)=2i\omega H_{2000}\left(\theta \right)-\Phi \left(\theta \right)f_{20000}^{\left(12\right)},\\ A\left(0\right)H_{1100}\left(\theta \right)=-\Phi \left(\theta \right)f_{1100}^{\left(12\right)},\\ A\left(0\right)H_{0200}\left(\theta \right)=-2i\omega H0200\left(\theta \right)-\varphi \left(\theta \right)f_{02000}^{\left(12\right)}\end{array}$ (4.22)

The initial condition of Equation (4.22) is also computed as

$\begin{array}{l}L{H}_{2000}=2i\omega H_{2000}\left(0\right)+f_{20000}^{\left(2\right)}-\Phi \left(0\right)f_{20000}^{\left(12\right)},\\ L{H}_{1100}=f_{11000}^{\left(2\right)}-\Phi \left(0\right)f_{11000}^{\left(12\right)},\\ L{H}_{0200}=-2i\omega H_{0200}\left(0\right)+f_{02000}^{\left(2\right)}-\Phi \left(0\right)f_{02000}^{\left(12\right)}\end{array}$ (4.23)

Based on the initial condition (4.23), one can derive the coefficients $H_{2000}\left(\theta \right)\mathrm{,}{H}_{1100}\left(\theta \right)\mathrm{,}{H}_{0200}\left(\theta \right)$ by the integral Equation (4.22) with respect to $\theta$. Therefore, we get the coefficient $k_{12}$ in formal norm (4.18) to be listed as

$\begin{array}{l}{k}_{12}=f_{21000,1}^{\left(1\right)}+\frac{f_{20000,1}^{\left(12\right)}{f}_{11000,1}^{\left(12\right)}}{i\omega }-\frac{{\left|f_{11000,1}^{\left(12\right)}\right|}^{\left(12\right)}}{i\omega }+\frac{2{\left|f_{02000,1}^{\left(12\right)}\right|}^2}{3i\omega }+f_{10100,y_0}^{\left(12\right)}{H}_{1100}\left(0\right)\\ +f_{10010}^{\left(12\right)}{H}_{1100}\left(-{\tau }_c\right)+f_{01100,y_0}^{\left(12\right)}{H}_{2000}\left(0\right)+f_{01010,y_{\tau }}^{\left(12\right)}{H}_{2000}\left(-{\tau }_c\right)\end{array}$ (4.24)

5. Conclusion

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 G₁/G₂ 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.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Murray, A. and Hunt, T. (1993) The Cell Cycle. Oxford University Press, Oxford.
[2]	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
[3]	Edelstein-Keshet, L. (1988) Mathematical Models in Biology. McGraW-Hill, Inc., New York.
[4]	Ukai, H. and Ueda, H.R. (2010) Systems Biology of Mammalian Circadian Clocks. Annual Review of Physiology, 72, 579-603. https://doi.org/10.1146/annurev-physiol-073109-130051
[5]	Qu,Z., MacLellan, W.R. and Weiss, J.N. (2003) Dynamics of the Cell Cycle: Checkpoints, Sizers, and Timers. Biophysical Journal, 85, 3600-3611. https://doi.org/10.1016/S0006-3495(03)74778-X
[6]	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. https://doi.org/10.1073/pnas.0308265100
[7]	Dubitzky, W. and Kurth, M.J. (2005) Mathematical Models of Cell Cycle Regulation. Bioinformatics, 6, 163-177. https://doi.org/10.1093/bib/6.2.163
[8]	Norel, R. and Agur, Z. (1991) A Model for the Adjustment of the Mitotic Clock by Cyclin and MPF Levels. Science, 251, 1076-1078. https://doi.org/10.1126/science.1825521
[9]	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. https://doi.org/10.1073/pnas.88.16.7328
[10]	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
[11]	Smith, H. (2011) An Introduction to Delay Differential Equations with Applications to the Life Sciences. Springer, New York.
[12]	Kuang, K. (1992) Delay Differential Equations with Application in Population Dynamics. Springer, New York.
[13]	Hale, J. (2003) Theory of Functional Differential Equations. World Publishing Corporation, Beijing.
[14]	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. https://doi.org/10.1142/S1793524510000829
[15]	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. https://doi.org/10.1142/S0218127409025080
[16]	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. https://doi.org/10.3934/dcdsb.2008.9.397
[17]	Xu, J. and Chung, K.W. (2003) Effects of Time Delayed Position Feedback on a Van Der Pol-Duffing Oscillator. Journal of Physics D, 180, 17-39. https://doi.org/10.1016/S0167-2789(03)00049-6
[18]	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. https://doi.org/10.1016/j.ijnonlinmec.2016.10.008