Hopf Bifurcation of a Gene-Protein Network Module with Reaction Diffusion and Delay Effects

Suqi Ma Ma

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

DOI: 10.4236/ijmnta.2021.103007                

Abstract

The infinite dimensional partial delay differential equation is set forth and delay difference state feedback control is considered to describe the cell cycle growth in eukaryotic cell cycles. Hopf bifurcation occurs as varying free parameters and time delay continuously and the multi-layer oscillation phenomena of the homogeneous steady state of a simple gene-protein network module is investigated. Normal form is derived based on normal formal analysis technique combined with center manifold theory, which is further to compute the bifurcating direction and the stability of bifurcation periodical solutions underlying Hopf bifurcation. Finally, the numerical simulation oscillation phenomena is in coincidence with the theoretical analysis results.

Keywords

Partial Functional Differential Equations, Hopf Bifurcation, Normal Form

1. Introduction

With the background knowledge of mathematical description of nonlinear dynamical model, people have endeavored to develop cell cycle growth model which manifests the mutual reaction among key components or its relationship with circumstance. Based on the fundamental tenets of cell biology that all the cells are derived from the preceding existing cells [1] [2], cells reproduce by duplicating their contents then divide into two daughter cells. During the cell cycle, two phases of inclusion should be emphasised: DNA synthesis stage (S phase), cell mitosis stage (M phase), and two “gaps” (G₀/G₁ phase) junction S phase and M phase. After chromosomes duplicate their contents, one mammaline cell divide into two daughter cells then the cell cycle completes its tasks to enter into the new cell cycle. The devastating thing is the existing models of the yeast cell cycle build up by quite detailed blocks yet [3] [4] [5]. People translate their acquired biology knowledge to some differential equations to describe the cell cycle growth in eukaryotic cell cycles [6] [7].

Yet as for now, the research works mainly focus on the functional sub-systems governed by gene-protein network modules since phase transition within it during cell cycle growth. Bifurcation analysis is helpful to identify the key components and their interaction relationship in the complex dynamical network. For example, refer to paper [8] [9], the authors put forward the double activator-inhibitor module to identify the positive feedback regulation in a mammalian gene-protein network to control G₁/S transition. In their paper, the bistability feature of steady states is observed and the authors highlight the system transition between stable steady state and its instability due to saddle-node bifurcation and also the transcritical bifurcation is considered.

As is well known, Michaelis-Menten rate law describes the activation feedback regulation function of gene-protein with the action of protein factors [10]. In addition, to drive the downstream events to generate gene protein production, enzymins reaction formed by binding itself to Cdk2 monomer. Time delay is incorporated into the phosphate groups while binding to target proteins to active protein phosphorylation process [11]. Therefore, the activation regulation with Cdk2 dimer dynamics is dominated by the hill function

$g\left(x\right)=\frac{s_0+x{\left(t-\tau \right)}^2}{{\theta }_1+x{\left(t-\tau \right)}^2}$

With the simple mathematical description, we put forth the following single gene-protein model in G₁/S phase transition,

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

wherein x denotes the gene-protein concentration, k₁ represents the degradation rate, and the second term describe the activation regulation function which acts on functional module. In addition, the last term denotes the nonlinear feedback control which depends on the concentration difference during junction gaps with the consideration of time delay, herein $2\tau$ means the total time of the sum of necessary time during G₀ and G₁ phase of the cell cycle.

The simulation work is finished by DDE-Biftool software which is applied to do dynamic analysis of delay differential equations with high technique. As shown in Figure 1(a), the bistability phenomena of the steady state of Equation (1.1) is observed. However, as varying free parameter k₁ and time delay, Hopf bifurcation further arise which change the stability property of the equilibrium solution. The bifurcating stable periodical solutions arising from the critical value of Hopf bifurcation point are continued as varying free parameter k₁ continuously. Without reaction diffusion, the periodical solution dies out when collide with the instable saddle and the homoclinic solution arise at the saddle with codimension

Figure 1. Form periodic solution to Homoclinic orbit arising in system (1.2) as varying free parameter k₁. (a) Stability and bifurcation of equilibrium solution and stability of the bifurcating periodic solution. (b) The time period of bifurcating periodic solutions are varied from finite period to infinity period. (c)The bifurcating homoclinic orbit when the periodic solution is collided with the unstable saddle. (d) The bifurcating homoclinic solution.

singularity 2. As shown in Figure 1(b), the time period of the bifurcating periodic solution is changed form finite time to infinity period. The bifurcating homoclinic orbit ${\Gamma }_i\left(t\right)$ is exhibited in Figure 1(c) and Figure 1(d) when collision phenomena happened which is satisfied by

$\underset{t\to \pm \infty }{lim}{\Gamma }_i\left(t\right)=S_i$

with $i=1,2$, and $S_1\mathrm{,}{S}_2$ are unstable saddles.

In this paper, the multi-layer oscillation phenomena are explored underlying Hopf bifurcation with diffusion effects. With the Neumman boundary condition, the bifurcating periodical waves are observed via varying free parameter underlying supercritical Hopf bifurcation. The continuation of periodical wave solution is also carried out in small parameter region.

Based on the fundamental theory of functional differential equation [12] [13] [14], people developed Lyapunov-Schmidt dimensional reduction scheme to compute the bifurcating direction of periodical solutions underlying Hopf bifurcation. Therefore, people applied center manifold theory in functional differential equation to compute the stability of bifurcating periodical solution [15] [16] [17] [18]. As for the partial functional differential equations, the well known center manifold theory has also been developed to further apply in normal form analysis near equilibrium solution and herein we adopt Teresa Faria’s method [19]. In Faria’s paper, based on theory of the autonomous functional differential equations, the analysis technique is addressed for calculating coefficients of normal form on center manifold.

The whole paper is organized as listed. In section 2, with homogenous Neumann conditions, the mathematical model of cell cycle growth model is described with reaction diffusion effects. In section 3, Hopf bifurcation is tracked as varying time delay and free parameter continuously. In section 4, based on the fundamental theory of partial functional differential equations, the normal form is computed with center manifold analytical technique, and finally the numerical simulation verifies the correctness of theoretical results.

2. The Mathematical Description with Diffusion Effects

With homogenous Neumann conditions, Equation (1.1) with free diffusion effect is modeled by

$\frac{\text{d}x}{\text{d}t}=D\frac{{\partial }^{2}x}{\partial a^2}-k_{1}x+k_2\frac{s_0+x{\left(t-\tau \right)}^2}{{\theta }_1+x{\left(t-\tau \right)}^2}+f\left(x\left(t-2\tau \right)-x\left(t\right)\right)$ (2.1)

with the definition $\Delta =\frac{{\partial }^2}{\partial a^2}$ is the Laplacian operator. Equation (2.1) is the

infinite dimensional functional partial differential equation defined on state space X and its definition domain satisfies $dom\left(\Delta \right)\subset X$. Without loss of generality, we assume X is the Hilbert space with inner product $\langle \varphi \mathrm{,}\psi \rangle ={\displaystyle {\int }_0^{2\pi }\varphi \left(x\right)\psi \left(x\right)\text{d}x}$. By setting $x\left(t+s\right)=x_t\left(s\right)$, the differntial operator of Equation (2.1) is defined on the infinite dimensional Banach space $C\left(\left[-2\tau \mathrm{,0}\right]\to X\right)$ with the super norm $\Vert \phi \Vert ={\sup }_{-2\tau \le \theta \le 0}\left|\phi \left(\theta \right)\right|$ for $\phi \in C$. The boundary and initial condition of Equation (2.1) is described as

$\begin{array}{l}\frac{\partial x}{\partial a}\left(\mathrm{0,}t\right)=\frac{\partial x}{\partial a}\left(2\pi \mathrm{,}t\right)=\mathrm{0,}a\in \left[\mathrm{0,2}\pi \right]\\ u_0=\varphi \left(\theta \right)\mathrm{,}\varphi \in C\left(\left[-2\tau \mathrm{,0}\right]\mathrm{,}X\right)\end{array}$ (2.2)

Hopf bifurcation occurs as varying free parameter and time delay, and the bifurcating periodical oscillating solutions are produced due to the instability phenomena. Hopf bifurcation occurs as the stability property of positive equilibrium solution first time changed at some critical value at some diffusion layer, Specially or not, we discuss Hopf bifurcation of system (2.1).

3. Hopf Bifurcation Analysis

Assume $E=x^{*}$ is the positive equilibrium solution of Equation (2.1) to satisfy

$-k_1{x}^{*}\left({\theta }_1+x^{*2}\right)+k_2\left(s_0+x^{*2}\right)=0$ (3.1)

Set $L\mathrm{<mo>\; :\; </mo>}C\left(\left[-2\tau \mathrm{,0}\right]\mathrm{,}R\right)\to X$ is the linear part of mobility. The nonlinear part $F{x}_t$ is the Taylor expansion beginning from quadratic term with $F\mathrm{<mo>\; :\; </mo>}C\left(\left[-2\tau \mathrm{,0}\right)\mathrm{,}R\right)\to X$. We define the complexification space

$X_C:=X\oplus iX=\left\{x_1+i{x}_2:x_1,x_2\in X\right\}$ (3.2)

The Taylor expansion of its truncation form of Equation (3.2) is written as

$\frac{\text{d}x}{\text{d}t}=D\left(\Delta x\right)+L{x}_t+F{x}_t$ (3.3)

It is verified that Equation (3.3) satisfies the following general condition:

(H1) $D\Delta$ generates a $C_0$ semigroup $\left\{T\left(t\right)\right\}$ on X with $\left|T\left(t\right)\right|\le M{\text{e}}^{wt}$ (for some $M\ge 1$ and $w\in R$ ) for all $t\ge 0$ ;

(H2) the eigenfunctions ${\left\{{\beta }_k\right\}}_{k=1}^{\infty }$ of $D\Delta$, generates orthonormal basis for X, and the corresponding eigenvalues ${\left\{{\mu }_k\right\}}_{k=1}^{\infty }$ satisfy ${\mu }_k\to -\infty$ ;

(H3) the subspaces $B_k\mathrm{<mo>\; :\; </mo>}\left\{{\langle v\mathrm{,}\beta \rangle }_k{\beta }_k\mathrm{<mo>\; |\; </mo>}v\in C\right\}$ of C satisfies $L\left(B_k\right)\subset span\left\{{\beta }_k\right\}$ ;

(H4) L can be extended to a bounded linear operator from BC to X wherein

$BC=\left\{\psi :C\left[-2\tau ,0\right)\to X|\psi \text{is}\text{continuous}\text{on}\left[-2\tau ,0\right),\exists {\lim }_{t\to 0^{-}}\psi \left(\theta \right)\in X\right\}$

with sup norm form.

Hopf Bifurcation

Based on the fundamental theory of partial functional differential equations as stated by [19], the linear differential operator (3.3) exists the unique solution which satisfies initial condition. And the generated strong continuous semigroup composed of solution operators has infinitesimal generator A:

$A\varphi =\stackrel{\dot{}}{\varphi }\left(\theta \right),dom\left(A\right)=\left\{\varphi \in C,\stackrel{\dot{}}{\varphi }\in C,\varphi \left(0\right)\in \left(dom\left(\Delta \right)\right),\stackrel{\dot{}}{\varphi }\left(0\right)=D\Delta \varphi \left(0\right)+L\varphi \right\}$

The operator A has only its point spectrum, with $\sigma \left(A\right)={\sigma }_P\left(A\right)=\left\{\lambda \in C|\Delta \left(\lambda \right)y=0,\text{with}y\in dom\left(A\right)-0\right\}$ and

$\Delta \left(\lambda \right)y=-\left({\mu }_m+\lambda \right)y+L\left({\text{e}}^{\lambda \cdot }\right)y$ (3.4)

It is well known that the eigenvalue problem

$-D{\varphi }^{″}=\mu \varphi ;a\in \left(0,2\pi \right),{\varphi }^{\prime }\left(0\right)={\varphi }^{\prime }\left(2\pi \right)=0$ (3.5)

has eigenvalues ${\mu }_m=D\frac{m^2}{4},m=0,1,2,3,\cdots$, with the corresponding eigenfunctions

${\varphi }_m=\cos \left(\frac{m}{2}x\right)$

Let $\varphi ={\displaystyle {\sum }_{i=0}^{\infty }}{a}_m{\varphi }_m$ be an eigenfunction of the eigenvalue problem (3.4), then we obtain a series of characteristic equation

$-{\mu }_m-\lambda +a+b{\text{e}}^{-\lambda \tau }+c{\text{e}}^{-2\lambda \tau }=0,m=0,1,2,3,\cdots$ (3.6)

Set $\lambda =i{\omega }_m\left({\omega }_m>0\right)$, then substitute it into Equation (3.6) and separate the real part from the imaginary part to get

$\begin{array}{l}\left(a+c-{\mu }_m\right)\cos \left({\omega }_m\tau \right)+{\omega }_m\sin \left({\omega }_m\tau \right)=-b,\\ \left(a-c-{\mu }_m\right)\sin \left({\omega }_m\tau \right)-{\omega }_m\cos \left({\omega }_m\tau \right)=0.\end{array}$ (3.7)

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

$\cos \left({\omega }_m\tau \right)=-\frac{\left(a-c-{\mu }_m\right)b}{{\left(a-{\mu }_m\right)}^2-c^2+{\omega }_m^2},\sin \left({\omega }_m\tau \right)=-\frac{b{\omega }_m}{{\left(a-{\mu }_m\right)}^2-c^2+{\omega }_m^2}$ (3.8)

Therefore, we have

${\left({\left(a-{\mu }_m\right)}^2-c^2+{\omega }_m^2\right)}^2-{\left(a-c-{\mu }_m\right)}^2{b}^2-b^2{\omega }_m^2=0$ (3.9)

Condition 1: ${\left(-a^2+2a{\mu }_m+\frac{1}{2}{b}^2+c^2-{\mu }_m^2\right)}^2\le \frac{1}{4}\left(-8a{b}^{2}c+b^4+8b^2{c}^2+8b^{2}c{\mu }_m\right)$, then Equation (3.8) has one pair of imaginary roots $\pm i{\omega }_m$ given that

${\omega }_m^2=-a^2+2a{\mu }_m+\frac{1}{2}{b}^2+c^2-{\mu }_m^2+\frac{1}{2}\sqrt{-8a{b}^{2}c+b^4+8b^2{c}^2+8b^{2}c{\mu }_m}$ (3.10)

Condition 2: ${\left(-a^2+2a{\mu }_m+\frac{1}{2}{b}^2+c^2-{\mu }_m^2\right)}^2>\frac{1}{4}\left(-8a{b}^{2}c+b^4+8b^2{c}^2+8b^{2}c{\mu }_m\right)$, then Equation (3.8) has two pair of imaginary roots $\pm i{\omega }_1\mathrm{,}\pm i{\omega }_2$ given that

${\omega }_{m1,2}^2=-a^2+2a{\mu }_m+\frac{1}{2}{b}^2+c^2-{\mu }_m^2\pm \frac{1}{2}\sqrt{-8a{b}^{2}c+b^4+8b^2{c}^2+8b^{2}c{\mu }_m}$ (3.11)

The critical time delay $\tau$ for Hopf bifurcation is

${\tau }_m=\frac{1}{{\omega }_m}\left(\arctan \frac{{\omega }_m}{a+c-{\mu }_m}+k\pi \right)$ (3.12)

for $k=0,1,2,\cdots$. With the aids of the above analysis, stability property for the positive equilibrium solution is plotted as shown in Figure 2(a), with the related parameter value $k_2=1$, $\theta =2$, $s=0.029$, $k_3=-0.6$, $\tau =1.2$. By varying free parameter $k_1$ and time delay $\tau$ continuously, Hopf lines are also pictured for $m=0,1,2,\cdots$. It exhibits that the critical value of free parameter for Hopf bifurcation satisfies $k_{1,0}>k_{1,1}>k_{1,2}\cdots$, hence Hopf bifurcation occurs at $k_1=k_{1,0}=0.4976$ while $m=0$.

Figure 2. Hopf bifurcation of the homogeneous equilibrium solution of system (3.3). (a) Stability property of equilibrium solution for $m=0,1,2$, wherein $\tau =1.2$ ; (b) Hopf lines on $\left(k_1\mathrm{,}\tau \right)$ plane for $m=0,1,2$.

4. Normal Form Computation

The periodical solution arise near Hopf point. Based on the known center manifold theory, people have applied dimensional reduction technique to analyze the bifurcating direction of periodical solution. Via the computation of the coefficients of norm form, we also explore the stability of periodical solutions. As is well known, the parameter perturbation scheme is useful in carrying out the computation of norm form coefficient to show the bifurcating direction of Hopf bifurcation.

Set $x=x-x^{*}$ with $x^{\mathrm{<mo>\; *\; </mo>}}$ is the unique positive equilibrium solution, we adopt the parameter perturbation method further to analyze Hopf bifurcation direction. With the assumption of $0<\epsilon \ll 1$, near Hopf point $\left(k_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)$, set $k_1=k_1^{*}+\epsilon k_e$, $\tau ={\tau }^{*}+\epsilon {\tau }_e$ then do dimensionless transformation $x\to \epsilon x$, one gets the abstract form of Equation (2.2) as

$\stackrel{\dot{}}{x}\left(t\right)=D\Delta x+L{x}_t+L_{\epsilon }{x}_t+F\left(x_t\right)$ (4.1)

wherein for $\varphi \in C$, there exists bounded variation function $\eta \mathrm{<mo>\; :\; </mo>}\left[-2\tau \mathrm{,0}\right]\to X$ which satisfy

$L\varphi ={\displaystyle {\int }_{-2{\tau }^{*}}^0}\text{d}\eta \left(\theta \right)\varphi \left(\theta \right)$ (4.2)

with

$\text{d}\eta \left(\theta \right)=a\delta \left(\theta \right)+b\delta \left(\theta +\tau \right)+c\delta \left(\theta +2\tau \right)$ (4.3)

and

$L_{\epsilon }\varphi =-\epsilon k_e\varphi +{\displaystyle {\int }_{-\tau }^0}\text{d}\eta \left(\theta \right)\varphi \left(\theta \right)-{\displaystyle {\int }_{-{\tau }^{*}}^0}\text{d}\eta \left(\theta \right)\varphi \left(\theta \right)$ (4.4)

In addition, we expand nonlinear part $F(\cdot )$ to be its Taylor form with 3^rd trunction as

$\begin{array}{c}F\left(\varphi \right)=-\frac{\left(2\left(-3x^{*2}+{\theta }_1\right)\right)k_2\left(s_0-{\theta }_1\right)}{{\left(x^{*2}+{\theta }_1\right)}^3}{\varphi }^2\left(-\tau \right)\\ -\frac{24k_2{x}^{*}\left(x^{*2}-{\theta }_1\right)\left(s_0-{\theta }_1\right)}{{\left(x^{*2}+{\theta }_1\right)}^4}{\varphi }^3\left(-\tau \right)\end{array}$ (4.5)

The linear version of Equation (4.1) is rewritten as

$\stackrel{\dot{}}{x}\left(t\right)=D\Delta x+L{x}_t$ (4.6)

And the generated strong continuous semigroup composed of solution operators has infinitesimal generator A:

$A\varphi =\stackrel{\dot{}}{\varphi }\left(\theta \right),dom\left(A\right)=\left\{\varphi \in C,\stackrel{\dot{}}{\varphi }\in C,\varphi \left(0\right)\in \left(dom\left(\Delta \right)\right),\stackrel{\dot{}}{\varphi }\left(0\right)=D\Delta \varphi \left(0\right)+L\varphi \right\}$ (4.7)

The corresponding adjoint operator $A^{\mathrm{<mo>\; *\; </mo>}}$ defined on the conjugate space $C^{\mathrm{<mo>\; *\; </mo>}}=C\left(\left[\mathrm{0,2}\tau \right]\mathrm{,}X\right)$ is written as

$\begin{array}{l}{A}^{*}\psi =\stackrel{\dot{}}{\psi }\left(s\right),\\ dom\left(A^{*}\right)=\left\{\psi \in C^{*},\stackrel{\dot{}}{\psi }\in C^{*},\psi \left(0\right)\in \left(dom\left(\Delta \right)\right),\stackrel{\dot{}}{\psi }\left(0\right)=-D\Delta \psi \left(0\right)-L^{*}\psi \right\}\end{array}$ (4.8)

with

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

Suppose $\Lambda =\left\{i{\omega }_k\mathrm{,}-i{\omega }_k\right\}$ is the set of eigenvalue with zero real parts for some $k\in N$, and other eigenvalues have negative real parts. Then set $B_k:=span{\beta }_k$, we define $L\left(\varphi \right){\beta }_k=L\left(\varphi {\beta }_k\right)$ and write Equation (4.6) as its equivalent form

$\stackrel{\dot{}}{z}\left(t\right)=-{\mu }_{k}z\left(t\right)+L_k{z}_t$ (4.10)

The adjoint bilinear form $\langle \varphi (\cdot )\mathrm{,}\psi (\cdot )\rangle$ on $C\times C^{\mathrm{<mo>\; *\; </mo>}}$, $\varphi \in C\mathrm{,}\psi \in C^{\mathrm{<mo>\; *\; </mo>}}$ is defined as

$\langle \psi ,\varphi \rangle =\psi \left(0\right)\varphi \left(0\right)+{\displaystyle {\int }_{-\tau }^0}\psi \left(\xi +\tau \right)b\varphi \left(\xi \right)\text{d}\xi +{\displaystyle {\int }_{-2\tau }^0}\psi \left(\xi +2\tau \right)c\varphi \left(\xi \right)\text{d}\xi$ (4.11)

Suppose $P_k$ is the eigensubspace corresponding to $\Lambda$, then the phase space C can be decomposed into

$\begin{array}{l}C=P_k\oplus Q_k,P_k=span\left\{{\Phi }_k\right\},P_k^{\mathrm{<mo>\; *\; </mo>}}=span\left\{{\Psi }_k\right\}\\ dim\left(P_k\right)=dim\left(P_k^{\mathrm{<mo>\; *\; </mo>}}\right)=\mathrm{2,}\langle {\Psi }_k\mathrm{,}{\Phi }_k\rangle =I\mathrm{,}{\stackrel{\dot{}}{\Phi }}_k={\Phi }_k{B}_k\end{array}$ (4.12)

with $B_k=\left(\begin{array}{cc}i{\omega }_k& 0\\ 0& -i{\omega }_k\end{array}\right)$. Define the projection operator $\Pi \mathrm{<mo>\; :\; </mo>}C\to P_k$ with

$P_k=Im\left(\Pi \right)\mathrm{,}C=Im\left(\Pi \right)\oplus Q_k$ (4.13)

herein, $Q_k$ is the complement subspace of $P_K$. For any $\varphi \in BC$, we can write $\varphi ={\varphi }_0+X_0\alpha$ with definition

$X_0=\{\begin{array}{ll}0,\hfill & -2\tau \le \theta <0,\hfill \\ 1,\hfill & \theta =0\hfill \end{array}$

For ${\varphi }_0\in C$, We also define $\Pi \mathrm{<mo>\; :\; </mo>}C\to P_k$

$\Pi \varphi ={\Phi }_k\langle {\Psi }_k\mathrm{,}\left(\varphi \mathrm{,}{\beta }_k\right)\rangle {\beta }_k$ (4.14)

and

$\Pi X_0\alpha ={\Phi }_k{\Psi }_k^{\text{T}}\left(0\right)X_0\left(\alpha \mathrm{,}{\beta }_k\right){\beta }_k$ (4.15)

Alike FDE reduction method, we want to enlarge the phase space in such a way that Equation (4.1) can be written as an anstract form of ODE on Banach space BC. For any $\varphi \in BC$, we write $\varphi ={\varphi }_0+X_0\alpha$ with definition $BC=C\times X$. BC is a Banach space with super norm $\left|\varphi \right|={\left|{\varphi }_0\right|}_C+{\left|\alpha \right|}_X$.

With the infinitesimal generator $A\varphi$ given in Equation (4.7), the extension of $A\mathrm{<mo>\; :\; </mo>}C\subset BC\to BC$ is written as

$A\varphi =\stackrel{\dot{}}{\varphi }+X_0\left(D\Delta \varphi \left(0\right)+L\varphi -\stackrel{\dot{}}{\varphi }\left(0\right)\right)$ (4.16)

with $\varphi \in C\mathrm{,}\stackrel{\dot{}}{\varphi }\in C\mathrm{,}dom\left(\Delta \right)\in X$, wherein $BC=\left\{\varphi |\varphi \left(\theta \right)\text{is}\text{continuous}\text{on}\left[-2\tau ,0\right),\exists {\lim }_{\theta \to 0^{-}}\varphi \left(\theta \right)\in X\right\}$. Similarly, the infinitesimal generator $A^{\mathrm{<mo>\; *\; </mo>}}$ defined by Equation (4.8) can be extended on $B{C}^{*}=C^{*}\times X$, but which is omitted here.

The projection leads to the decomposition of the extended phase space as

$BC=P\oplus Ker\left(\Pi \right)$ (4.17)

with the property $Q⊊Ker\left(\Pi \right)$, and linear operator A commutes with operator $\Pi$.

Set $x_t=u\left(t\right)$, Equation (4.1) can be written as its abstract ODE form on the extended phase space BC,

$\frac{\text{d}}{\text{d}t}u\left(t\right)=Au\left(t\right)+A_{\epsilon }u\left(t\right)+R\left(u\left(t\right)\right)$ (4.18)

with

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

Set $\tilde{z}\left(t\right)=\left({\Psi }_k\mathrm{,}\langle \varphi \mathrm{,}{\beta }_k\rangle \right)$ and $z\left(t\right)={\Phi }_k\tilde{z}\left(t\right){\beta }_k+y_t{\beta }_k$, then the linear part is transformed into the following form,

$\begin{array}{c}{\tilde{z}}^{\prime }\left(t\right)=B_k\tilde{z}\left(t\right)+{\Psi }_k^{\text{T}}\Pi X_0{A}_{\epsilon }\left({\Phi }_k\tilde{z}\left(t\right){\beta }_k+y\left(t\right){\beta }_k\right)\\ =B_k\tilde{z}\left(t\right)+{\Psi }_k^{\text{T}}\left(0\right)\langle L_{\epsilon }\left({\Phi }_k\tilde{z}\left(t\right){\beta }_k\right),{\beta }_k\rangle \\ =B_k\tilde{z}\left(t\right)+{\Psi }_k^T\left(0\right)\left(-\epsilon k_e{\Phi }_k\left(0\right)-\epsilon b{B}_k{\Phi }_k\left(-{\tau }^{*}\right){\tau }_e\right)\tilde{z}\left(t\right)\end{array}$ (4.19)

Further, considering the nonlinear part $F(\cdot )$, the dimensional reduction system of Equation (4.1) is written as

$\begin{array}{l}{\tilde{z}}^{\prime }\left(t\right)=B_k\tilde{z}\left(t\right)+{\Psi }_k^{\text{T}}\Pi X_0\langle R\left({\Phi }_k\tilde{z}\left(t\right){\beta }_k+y_t{\beta }_k\right),{\beta }_k\rangle \\ =B_k\tilde{z}\left(t\right)+{\Psi }_k^{\text{T}}\left(0\right)\langle F\left({\Phi }_k\tilde{z}\left(t\right){\beta }_k+y_t{\beta }_k\right),{\beta }_k\rangle \\ y^{\prime }\left(t\right)=Ay\left(t\right)+\left(I-\Pi \right)X_0\langle R\left({\Phi }_k\tilde{z}\left(t\right){\beta }_k+y_t{\beta }_k\right),{\beta }_k\rangle \\ =Ay\left(t\right)+\{\begin{array}{l}-{\Phi }_k\left(\theta \right){\Psi }_k^{\text{T}}\left(0\right)\langle F\left({\Phi }_k\tilde{z}\left(t\right){\beta }_k+y_t{\beta }_k\right),{\beta }_k\rangle \hfill \\ \langle F\left({\Phi }_k\tilde{z}\left(t\right){\beta }_k+y_t{\beta }_k\right),{\beta }_k\rangle -{\Phi }_k\left(0\right){\Psi }_k^{\text{T}}\left(0\right)\langle F\left({\Phi }_k\tilde{z}\left(t\right){\beta }_k+y_t{\beta }_k\right),{\beta }_k\rangle \hfill \end{array}\end{array}$ (4.20)

Suppose

$\begin{array}{c}\hat{f}={\Psi }_k^{\text{T}}\left(0\right)\langle F\left({\Phi }_k\tilde{z}\left(t\right){\beta }_k+y{\beta }_k\right),{\beta }_k\rangle \\ =f_{200}^{\left(1\right)}{z}_1^2+f_{110}^{\left(1\right)}{z}_1{z}_2+f_{020}^{\left(1\right)}{z}_2^2+f_{101}^{\left(1\right)}{z}_{1}y\left(-\tau \right)+f_{011}^{\left(1\right)}{z}_{2}y\left(-\tau \right)+f_{210}^{\left(1\right)}{z}_1^2{z}_2+\cdots \end{array}$ (4.21)

with

$\begin{array}{l}{f}_{200}^{\left(1\right)}=\frac{1}{\pi \sqrt{\pi }}{s}_1{\Psi }_k^{\text{T}}\left(0\right){\Phi }_{k1}\left(-{\tau }^{*}\right)\left({\left(-1\right)}^m+{\left(-1\right)}^{3m}-2\right),\\ f_{110}^{\left(1\right)}=\frac{2}{\pi \sqrt{\pi }}{s}_1{\Psi }_k^{\text{T}}\left(0\right){\Phi }_{k1}\left(-{\tau }^{*}\right){\Phi }_{k2}\left(-{\tau }^{*}\right)\left({\left(-1\right)}^m+{\left(-1\right)}^{3m}-2\right),\\ f_{020}^{\left(1\right)}=\frac{1}{\pi \sqrt{\pi }}{s}_1{\Psi }_k^{\text{T}}\left(0\right){\Phi }_{k2}\left(-{\tau }^{*}\right)\left({\left(-1\right)}^m+{\left(-1\right)}^{3m}-2\right),\\ f_{101}^{\left(1\right)}=\frac{2}{\pi \sqrt{\pi }}{s}_1{\Psi }_k^{\text{T}}\left(0\right){\Phi }_{k1}\left(-{\tau }^{*}\right)\left({\left(-1\right)}^m+{\left(-1\right)}^{3m}-2\right),\\ f_{011}^{\left(1\right)}=\frac{2}{\pi \sqrt{\pi }}{s}_1{\Psi }_k^{\text{T}}\left(0\right){\Phi }_{k2}\left(-{\tau }^{*}\right)\left({\left(-1\right)}^m+{\left(-1\right)}^{3m}-2\right),\\ f_{210}^{\left(1\right)}=\frac{9}{4\pi }{s}_2{\Phi }_{k1}^2\left(-{\tau }^{*}\right){\Phi }_{k2}\left(-{\tau }^{*}\right)\end{array}$

Furthermore, we set

$y\left(t\right)=H_{20}\left(\theta \right)z_1^2+H_{11}\left(\theta \right)z_1{z}_2+H_{02}\left(\theta \right)z_2^2+\cdots$ (4.22)

Substitute it into Equation (4.17), one obtains

$\begin{array}{l}A{H}_{20}\left(\theta \right)=2i\omega H_{20}\left(\theta \right)-{\Phi }_k\left(\theta \right)f_{200},\\ A{H}_{11}\left(\theta \right)=-{\Phi }_k\left(\theta \right)f_{110},\\ A{H}_{02}\left(\theta \right)=-2i\omega H_{02}\left(\theta \right)-{\Phi }_k\left(\theta \right)f_{020},\end{array}$ (4.23)

with the initial value condition

$\begin{array}{l}-{\mu }_{k}D{H}_{20}+L{H}_{20}=2i\omega H_{20}\left(0\right)-{\Phi }_k\left(0\right)f_{200}^1+f_{200}^{\left(2\right)},\\ -{\mu }_{k}D{H}_{11}+L{H}_{11}={\Phi }_k\left(0\right)f_{110}^{\left(1\right)}+q_{110}^{\left(2\right)},\\ -{\mu }_{k}D{H}_{02}+L{H}_{02}=-2i\omega H_{02}\left(0\right)-{\Phi }_k\left(0\right)f_{020}^{\left(1\right)}+f_{020}^{\left(2\right)},\end{array}$ (4.24)

with

$\begin{array}{l}{f}_{200}^{\left(2\right)}=\frac{1}{\pi \sqrt{\pi }}{s}_1{\Phi }_{k1}\left(-{\tau }^{*}\right)\left({\left(-1\right)}^m+{\left(-1\right)}^{3m}-2\right),\\ f_{110}^{\left(2\right)}=\frac{2}{\pi \sqrt{\pi }}{s}_1{\Phi }_{k1}\left(-{\tau }^{*}\right){\Phi }_{k2}\left(-{\tau }^{*}\right)\left({\left(-1\right)}^m+{\left(-1\right)}^{3m}-2\right),\\ f_{020}^{\left(2\right)}=\frac{1}{\pi \sqrt{\pi }}{s}_1{\Phi }_{k2}\left(-{\tau }^{*}\right)\left({\left(-1\right)}^m+{\left(-1\right)}^{3m}-2\right)\end{array}$

By the near identity transformation, we obtain that

${\tilde{z}}^{\prime }\left(t\right)=B_k\tilde{z}\left(t\right)+\epsilon \tilde{b}\tilde{z}\left(t\right)+\tilde{d}{\tilde{z}}^2\overline{\stackrel{˜}{z}}$ (4.25)

with

$\begin{array}{l}\tilde{b}={\Psi }_k\left(0\right)\left(-k_e{\Phi }_k\left(0\right)-b{B}_k{\Phi }_k\left(-{\tau }^{*}\right){\tau }_e\right)\\ \tilde{d}=f_{210}^{\left(1\right)}+\frac{f_{200}^{\left(1\right)}{f}_{110}^{\left(1\right)}}{i{\omega }_k}-\frac{f_{110}^{\left(1\right)}}{i{\omega }_k}+\frac{2f_{020}^{\left(1\right)}}{3i\omega }+f_{101}^{\left(1\right)}{H}_{11}\left(-{\tau }^{*}\right)+f_{011}^{\left(1\right)}{H}_{20}\left(-{\tau }^{*}\right)\end{array}$ (4.26)

Note that herein, we suppose that the multiplication of the vector $v_1=\left(\begin{array}{c}{a}_1\\ b_1\end{array}\right)$ and $v_2=\left(\begin{array}{c}{a}_2\\ b_2\end{array}\right)$ means the multiplication between row elements, that is $v_1{v}_2=\left(\begin{array}{c}{a}_1{a}_2\\ b_1{b}_2\end{array}\right)$.

Based on the above analysis, we have the following theorem,

Theorem 3.1. The norm form of Equation (3.1) near Hopf point $\left(k_1^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{\tau }^{\mathrm{<mo>\; *\; </mo>}}\right)$ can be written as

$z\left(t\right)=i\omega z\left(t\right)+\epsilon {\tilde{b}}_{11}z\left(t\right)+{\tilde{d}}_1{z}^2\bar{z}$ (4.27)

Hence, the periodical solutions with small amplitude arise underlying Hopf bifurcation if $\mu =\frac{Re\left({\tilde{d}}_1\right)}{Re\left({\tilde{b}}_{11}\right)}<0$, and the bifurcating solution is stable if $Re\left({\tilde{d}}_1\right)<0$, and unstable on reverse.

For example, it is calculated that Hopf bifurcation occurs at Hopf point $E_1=\left(k_1,\tau \right)=\left(0.258,1\right)$ and $E_2=\left(k_1,\tau \right)=\left(0.2168,1\right)$ while $k=0$. Respectively, the stable bifurcating periodical solution arise at Hopf point $E_1\mathrm{,}{E}_2$ which is supercritical. Near $E_1$, periodical oscillating solutions are computed respectively with $D=0.00001$ and $D=0.00002$. As shown in Figure 3(b), the periodical solution with maximal and minimal amplitudes are simulated by chosing $n=100$ to denote diffusion layer, time step $h=0.01$. In Figure 3(a), the bistable

Figure 3. The equilibrium solutions and the amplitude of bifurcating periodical solutions at Hopf point E₁. (a) The equilibrium solution continuously with free parameter k₁ and Hopf bifurcation occurs at $E_1=\left(k_1,\tau \right)=\left(0.258,1\right)$ and $E_2=\left(k_1,\tau \right)=\left(0.2168,1\right)$. (b) The maximal and minimal amplitudes of the bifurcating periodical solutions with free parameter k₁ varying continuously, while the diffusion coefficient is $D=0.00001$ and $D=0.00002$ respectively.

Figure 4. The bifurcating periodical solution of Hopf singularity. (a) The equilibrium is asymptotically stable; (b) The observed periodical solution induced by Hopf bifurcation; (c) The time series solutions with diffusion coefficient $D=0.00001$ and $D=0.00002$ ; (d) The phase portraits with diffusion coefficient $D=0.00001$ and $D=0.00002$.

coexistence of equilibrium solutions are observed and Hopf bifurcation occurs at $E_1\mathrm{,}{E}_2$. In Figure 4, the observed oscillating periodical solutions are induced due to the instability transition of equilibrium solution. As shown in Figure 4(a) and Figure 4(b), the equilibrium solution is asymptotically stable with $k_1=0.2758$, however Hopf bifurcation at $E_1$ leads to the stable oscillating periodical solution bifurcation. The periodical solution is observed on some diffusion layer, and the time series solution and the corresponding solution are shown in Figure 4(c) and Figure 4(d) respectively. Near $E_2$, the periodical solution is also observed as shown in Figure 5(b) which is bifurcated from the corresponding equilibrium solution. The equilibrium solution is asymptotically stable as observed in Figure 5(b). The interesting phenomena of continuous periodical

Figure 5. The bifurcating periodical solution arise from Hopf bifurcation point. (a) The equilibrium is asymptotically stable; (b) The observed periodical solution induced by Hopf bifurcation; (c) The time series solutions of two different diffusion layers; (d) The phase portraits of two different diffusion layers.

oscillation with different initial phase is observed in Figure 5(a). By chosen $n=10$ and $n=50$, the phase portraits and time series solutions are observed on two diffusion layers, as shown in Figure 5(c) and Figure 5(d).

5. Conclusion

The partial delay differential equation of gene reaction protein equation was set forth. The stability dynamics and Hopf bifurcation was analyzed underlying the feedback control of state difference between present state and its past time state. Without diffusion effects, using DDE-Biftool software, the bifurcating homoclinc

Figure 6. The continuation of homoclinic orbit as varying free parameter and time delay. (a) The bifurcating homoclinc solutions due to collison phenomena of continuation periodic solutions form Hopf point E₁ with the unstable saddle; (b) The bifurcating homoclinc solutions by collision of continuation of periodical solution bifurcated from Hopf point E₂ with the saddle.

solution was derived as the time period tends to infinity with continuation of periodic solution as varying free parameter. The continuation of homoclinc orbit becomes a possible job with application of DDE-Biftool software, as shown in Figure 6(a) and Figure 6(b). With diffusion effects, Hopf bifurcation phenomena were further analyzed and the multi-layer periodical oscillation phenomena were discovered. Combined with center manifold technique, the bifurcating direction of periodical solution was analyzed with norm form analysis method.

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]	Tyson, J.J., Novak, B., Odell, G.M., et al. (1996) Chemical Kinetic Theory: Understanding Cell-Cycle Regulation. Trends in Biochemical Sciences, 21, 89-96. https://doi.org/10.1016/S0968-0004(96)10011-6
[3]	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
[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. Briefings in 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, 3733-3751. https://doi.org/10.1142/S0218127409025080
[16]	Ma, S.Q., Feng, Z.S. and Lu, Q.S. (2008) A Two Parameter Criteria for Delay Differential Equations. Discrete & Continuous Dynamical Systems-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. Physica D: Nonlinear Phenomena, 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
[19]	Teresa, F. (2000) Normal Forms and Hopf Bifurcation for Partial Diffeential Equations with Delays. Transactions of the American Mathematical Society, 352, 2217-2238. https://doi.org/10.1090/S0002-9947-00-02280-7