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

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$.

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

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

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

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.