# Stability and Hopf Bifurcation Analysis of a Gene Expression Model with Diffusion and Time Delay

Authors: Yahong Peng; Tonghua Zhang
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We consider a model for gene expression with one or two time delays and diffusion. The local stability and delay-induced Hopf bifurcation are investigated. We also derive the formulas determining the direction and the stability of Hopf bifurcations by calculating the normal form on the center manifold.


---

## Body

## 1. Introduction

The study on dynamics of a biological model is one of the dominant subjects in mathematical biology due to its universal existence and importance. In this paper, we consider a mathematical model of intracellular regulatory system which began with the work of Goodwin [1]. Since then, many researchers developed this work [2–5]. However, it is Monk who developed the first mathematical model for the Hes1 system and validated it with biological data [6]. In order to describe the intracellular process more precisely, he introduced time delays to account for the processes of transcription and translation. If we denote by M(t) and P(t) the concentrations of Hes1 mRNA and Hes1 protein, respectively, the basic reaction kinetics for this system can be expressed in the form of
(1)dM(t)dt=f(P(t-τm))-cM(t),dP(t)dt=aM(t-τp)-bP(t),
where f(P)=α/(1+(P/P0)n). The parameter a is the rate at which Hes1 protein is produced from Hes1 mRNA and b and c are the decay rates of Hes1 protein and Hes1 mRNA, respectively. f(P) is the rate of production of new mRNA molecules, with α and P0 as constants, to represent the rate of transcript initiation in the absence of Hes1 protein and the reference concentration of protein, respectively, and n is the Hill coefficient. τm and τp represent the transcriptional and translational time delays. The units of the parameters are as follows. a is measured in protein molecules per mRNA molecule per minute; b and c are measured in molecules per minute; α is measured in mRNA molecules per diploid cell per minute; P0 is measured in molecules and time delays τm and τp are measured in minute.In order to reduce the number of parameters, [7] introduced transformations
(2)m=Mα,p=Pαa,p0=P0αa,
under which system (1) takes the following form:
(3)dm(t)dt=11+(p(t-τm)/p0)n-cm(t),dp(t)dt=m(t-τp)-bp(t).
Recently, Zhang et al. [8] investigated the stability and Hopf bifurcation of the equilibrium of system (3). By using the method of multiple time scales, they also obtained the normal form on the center manifold of delay differential equations (3).The diffusion process comes naturally in biology as pointed out by Murray in [9]. Experimental data also suggests that many pathways exhibit oscillation in concentrations of substance involved, both temporally and spatially. For example, a type diffusion process appears in a process of assemblage of particles, which is due to the particles spread out as a result of an irregular individual particle's motion. In these cases, a partial differential equation or system of partial differential equations will be employed to describe the processes. Nowadays, models involving delays and also spatial diffusion are increasingly applied to the study due to more appropriate biological justification. Reference [10] gave plausible biological explanations for the delays appearing in the model (3). In what follows, spatial diffusion will be applied to the model (3) to gain new information about the precise spatiotemporal dynamics of mRNA and proteins. In fact, Sturrock et al. [11] derived a system of partial differential equations to capture the evolution in space and time of the variables in the Hes1 and p53-Mdm2 systems. In this paper, we extend previous mathematical model (3) into its spatial version by considering spatial interactions within the cell explicitly (of course, the clear biological explanation will be given behind the model (4) of promotion); namely, we consider the delayed reaction-diffusion system with the following initial and boundary conditions:
(4)∂m(x,t)∂t=d1∂2m(x,t)∂x2+11+(p(x,t-τm)/p0)n-cm(x,t),t>0,x∈(0,π),∂p(x,t)∂t=d2∂2p(x,t)∂x2+m(x,t-τp)-bp(x,t),00000000000000000000000000t>0,x∈(0,π),∂m(x,t)∂x=∂p(x,t)∂x=0,x=0,π,t≥0,m(x,t)=ϕ(x,t)≥0,p(x,t)=ψ(x,t)≥0,0000000000000000000(x,t)∈[0,π]×[-τ,0],
where d1 and d2 are the diffusion coefficients for Hes1 mRNA and Hes1 protein, respectively, with unit such as cm/min. The initial function ϕ(t,x),ψ(t,x)∈𝒞:=C([-τ,0],L2([0,π])). The imposed Neumann boundary condition here implies that mRNA and protein are not exported across the nuclear membrane or the cell membrane.Obviously, (4) is a system of reaction-diffusion equations modeling the spatiotemporal evolution of the Hes1 system. The same reaction kinetics from the ODE model (3) are retained but are now also coupled with diffusion to model explicitly protein and mRNA transport in a cell. That is to say, molecules move from the nucleus to the cytoplasm and from cytoplasm to nucleus across the nuclear membrane. Here, we use a system of PDEs (4) to reflect the reality that the mRNA is transcribed from DNA exclusively in the nucleus and that protein is translated from mRNA exclusively in the cytoplasm. The main advantage of using systems of PDEs (4) to model intracellular reactions is that the PDEs enable spatial effects to be examined explicitly. The main object of this paper is to investigate the effect of the delay and diffusion on the dynamics of system (4). In addition, in order to determine the direction and the stability of Hopf bifurcations, we use the normal form procedure for functional differential equations (FDEs) with diffusion due to Faria [12, 13].This paper is organized as follows. In Section2, stability of positive equilibrium and existence of Hopf bifurcation are studied using τ as a parameter. In Section 3, the effect of diffusion on the Hopf bifurcation will be investigated. Using the normal form technique for partial functional differential equations, the formulas for determining the direction and stability of Hopf bifurcation are presented in Section 4. Finally, in Section 5, we will illustrate the theoretical results by numerical simulations along with some discussion.


## 2. Stability of Positive Equilibrium and Existence of Hopf Bifurcation

It is easy to verify that system (4) has a unique positive equilibrium E*(m*,p*) determined by
(5)11+(p*/p0)n-cm*=0,m*-bp*=0.The object here is then to relate the dynamics of (3) and (4) in the neighborhood of E*, at the first critical point of the parameter τ. To this end, we let m-=m-m*, p-=p-p*. With the help of (5), after dropping the bars for simplicity of notation, (4) is transformed into the following system:
(6)∂m(x,t)∂t=d1∂2m(x,t)∂x2+11+((p(x,t-τm)+p*)/p0)n-c(m(x,t)+m*),∂p(x,t)∂t=d2∂2p(x,t)∂x2+m(x,t-τp)-bp(x,t)
with the origin as its equilibrium.Let(7)f(1)(m,p)=11+((p(x,t-τm)+p*)/p0)n-c(m(x,t)+m*),f(2)(m,p)=m(x,t-τp)-bp(x,t).
Define fij(1)(i+j≥1) and fij(2)(i+j≥1) as follows:
(8)fij(1)=∂i+jf(1)(0,0)∂mi∂pj,fij(2)=∂i+jf(2)(0,0)∂mi∂pj,
where i and j are the nonnegative integers. Then system (6) can be rewritten as
(9)∂m(x,t)∂t=d1∂2m(x,t)∂x2+a11m(x,t)+a12p(x,t-τm)+∑i+j≥2‍1i!j!fij(1)mi(x,t)pj(x,t-τm),∂p(x,t)∂t=d2∂2p(x,t)∂x2+a21m(x,t-τp)+a22p(x,t)+∑i+j≥2‍1i!j!fij(2)mi(x,t-τp)pj(x,t),
where
(10)a11=f10(1)=∂f(1)∂m(0,0)=-c,a12=f01(1)=∂f(1)∂p(0,0)=-nc2b2(p*)n+1p0n,a21=f10(2)=∂f(2)∂m(0,0)=1,a22=f01(2)=∂f(2)∂p(0,0)=-b.
For simplification of notation, we use m(t) for m(·,t) and p(t) for p(·,t) and (m(t),p(t))=(m(·,t),p(·,t)) is in a suitable Hilbert space X(11)X={(m,p):m,p∈W2,2(0,π),∂m∂x=∂p∂x=0atx=0,π}.
By setting U(t)=(m(t),p(t))∈X, we further write (9) as an abstract equation in 𝒞:=C([-τ,0],X):
(12)dU(t)dt=dΔU(t)+L(Ut)+F(Ut),
where dΔ=(d1Δ,d2Δ) and L:𝒞→R2, and F:𝒞×R→R2 are given by
(13)L(φ)=(a11φ1(0)+a12φ2(-τm)a21φ1(-τp)+a22φ2(0)),F(φ)=(∑i+j≥2‍1i!j!fij(1)φ1i(0)φ2j(-τm)∑i+j≥2‍1i!j!fij(2)φ1i(-τp)φ2j(0)),
for φ=(φ1,φ2)∈𝒞. Obviously, L is a linear operator. The linearization of (12) is
(14)dU(t)dt=dΔU(t)+L(Ut).
It has a characteristic equation given by
(15)Δk(λ,τ)=λ2+pkλ+rk+Ke-2λτ=0,k=0,1,2,…,
where pk=(d1+d2)k2+(b+c), rk=(d1k2+c)(d2k2+b), K=nc2b2(p*)n+1/p0n, and 2τ=τm+τp is the total time delay. When there is no diffusion effect, namely, d1=d2=0, (15) can be written as
(16)λ2+(b+c)λ+bc+Ke-2λτ=0,
which is equivalent to (12) of Zhang et al.’s work in [8], and the related stability and Hopf bifurcation have been investigated.In what follows, we will analyze the effect of diffusion terms by the distribution of the roots of (15) with d1>0, d2>0. We first consider the case when the delay is zero. For (15), if τ=0, then we have
(17)λ2+pkλ+rk+K=0,k=0,1,2.
Since pk>0, rk>0 for any k∈N and K>0, it is easy to verify that (17) has a pair of roots with negative real parts. And, for τ>0, we have the following lemma.Lemma 1.
Assume that(H)(d1+c)(d2+b)≥K
holds. Then all the roots of the characteristic equation (15) have negative real part for τ>0.Proof.
If the conclusion is not true, namely, (15) admits at least one root λ=μ+iω with μ≥0, then we obtain
(18)(μ+iω)2+pk(μ+iω)+rk+Ke-2τ(μ+iω)=0.
Separating the real and imaginary parts yields
(19)μ2-ω2+pkμ+rk+Ke-μ·2τcos⁡(ω·2τ)=0,2μω+pkω-Ke-μ·2τsin(ω·2τ)=0.
It implies that
(20)(μ2-ω2+pkμ+rk)2+(2μω+pkω)2=K2e-2μ·2τ;
namely,
(21)μ4+ω4+pk2μ2+rk2+2μ2ω2+2pkμ3+2μ2ω2+(pk2-2rk)ω2+2pkμrk=K2e-2μ·2τ.
Since pk2-2rk=k4(d12+d22)+2(bd2+cd1)k2+b2+c2>0, we can easily verify
(22)rk2<K2e-2μ·2τ<K2.
That is, rk<K. Notice that (d1+c)(d2+b)≤rk, for k≥1, k∈N, gives
(23)(d1+c)(d2+b)<K.
It is a contradiction. Thus the conclusion follows.Notice that whenτ=0, all roots of (15) have negative real part and the roots of (15) continuously depend on the parameter τ, and we can summarize our conclusion as follows.Theorem 2.
Assume that(H) holds. Then the equilibrium point E* of system (4) is asymptotically stable for τ≥0.


## 3. Effects of Diffusion on the Hopf Bifurcation

Assume thatiω(ω>0) is a purely imaginary root of (15); then we have
(24)-ω2+rk+Kcos⁡(ω·2τ)=0,pkω-Ksin(ω·2τ)=0,
which implies that
(25)ω4+Pkω2+Rk=0,k=0,1,2,…,
where
(26)Pk=pk2-2rk=(d12+d22)k4+2(bd2+cd1)k2+b2+c2>0,Rk=rk2-K2=(rk+K)(rk-K)=[(d1k2+c)(d2k2+c)+K]×[(d1k2+c)(d2k2+c)-K].
Equation (25) implies that Rk should be negative for some k∈N. It is equivalent to the fact that
(27)R~k=(d1k2+c)(d2k2+c)-K
should be less than 0 as it is easy to see that Rk has the same sign as that of R~k. Rewrite R~k into the following form:
(28)R~k=d1d2k4+(d1b+d2c)k2+bc-K.
It is obviously a quadratic polynomial in terms of k2. Equation (28) implies that there exists a k1∈N such that R~k1<0 if and only if bc-K<0. Furthermore, we have
(29)R~k<0for0≤k≤k1,R~k>0fork>k1,k∈ℕ.
From (29), we obtain that, for each k∈{0,1,…,k1}, (25) has only one positive real root ωk, which is given by
(30)ωk=22-Pk+Pk2-4Rk.We now can make the following conclusion.Lemma 3.
Assume that(H) is not true. Then (15) has a pair of purely imaginary roots ±iωk for each k∈{0,1,…,k1} and has no purely imaginary roots for k1<k∈N, where k1 and ωk are defined as above.In the rest of this section, we will discuss the case ofk<k1.From (24), we have
(31)sin(ω·2τ)=pkωK,cos⁡(ω·2τ)=ω2-rkK.
Then, for k∈{0,1,…,k1}, define
(32)τjk=1ωk(arccosωk2-rkK+2jπ),j=0,1,2,….
In the following, we will order the sequence of τjk depending on the diffusion coefficients d1 and d2 for k∈{0,1,…,k1}.Lemma 4.
Ifd1=d2 and (H-) holds, then
(33)τ00=min⁡{τjk}k∈{0,1,…,k1},j=0,1,2,….Proof.
Ifd1=d2, from (30), we have
(34)ωk2=12[(b-c)2[2dk2+(b+c)]2+4K2-(b+c)2-2(b+c)dk2-2d2k4+2bci+(b-c)2[2dk2+(b+c)]2+4K2],ωk2-rk=12[+(b-c)2[2dk2+(b+c)]2+4K2-(b+c)2-4(b+c)dk2-4d2k4i+(b-c)2[2dk2+(b+c)]2+4K2].
Let x=(b-c)2[2dk2+(b+c)]2+4K2; it is easy to verify
(35)x>(b-c)2,(36)ωk2=12[x-x2-4K22(b-c)2+(b-c)22],ωk2-rk=12[x-x2-4K22(b-c)2].
Thus, according to (32), we obtain
(37)τjk=τ(x)=arccos[(x-((x2-4K2)/(b-c)2)/2K)+2jπ](2/2)[x-((x2-4K2)/2(b-c)2)+((b-c)2/2)]1/2.
In addition,
(38)d(x-((x2-4K2)/2(b-c)2))dx=1-x(b-c)2.
From (35), we know that x-((x2-4K2)/2(b-c)2) is decreasing with respect to x. Then we obtain τjk+1>τjk if d1=d2. Clearly, τj+1k>τjk, so we have
(39)τ00=min⁡{τjk}k∈{0,1,2,…},j=0,1,2,….According to Lemma4 and the continuous dependence of τjk on d1 and d2, we summarize the following lemma.Lemma 5.
For anyd>0, there exists an ϵ(d)>0 such that, for any d1,d2∈(d-ϵ,d+ϵ) satisfying (H-),τ00=min⁡{τjk}k∈{0,1,…,k1},j=0,1,2,….Letλk(τ)=μk(τ)+iωk(τ) be the roots of (15) near 2τ=τjk satisfying μk(τjk)=0, ωk(τjk)=ωk. By using the method in [14, 15], we can prove the following transversality condition.Lemma 6.
Ifd1 and d2 satisfy the condition in Lemma 5, then, for k∈{0,1,…,k1} and j∈ℕ0, dRe(λ)/dτ|2τ=τjk>0.Proof.
Differentiating equation (15) with respect to τ, we obtain
(40)(dλdτ)-1=(2λ+pk)e2λτ2λK-τλ.
From (24), we have
(41)Re(dλdτ|2τ=τjk)-1=2ωkcos⁡(ωkτjk)+pksin(ωkτjk)2kωk=2ωk·((ωk2-rk)/K)+pk·(pkωk/K)2kωk=2ωk2+(pk2-2rk)2K2>0
since pk2-2rk>0.Combining the above analysis and the qualitative theory of partial functional differential equations in [16], we have the following results on the stability of equilibrium E* of system (4) and existence of Hopf bifurcation near E*.Theorem 7.
Assume thatd1, d2 satisfy the condition in Lemma 5. Then, for system (4),(i)
the positive equilibriumE* is asymptotically stable for τ∈[0,τ00) and unstable for τ∈(τ00,+∞);(ii)
it undergoes Hopf bifurcations near the positive equilibriumE* at 2τ=τjk for k∈{0,1,…,k1} and j∈ℕ0.Remark 8.
Ifd1=d2=0, Theorem 2 and Theorem 7 are the conclusion of Theorem 1 in [8]. That is, assuming either condition (H) or (H-), Theorem 2 and Theorem 7 show that the local stability of E* for 0≤τ<τ00 is the same for system (3) and system (4). Here we can know the effect of the diffusion coefficients d1 and d2, or, in order to more clearly understand the effect of diffusion, we take b=c=0.03, p0=0.4, n=2, and then we get p*=5.6134, K≈17.91×10-4 by calculation. In the absence of diffusion, b×c=9×10-4<K, it is known that the equilibrium E* of system (3) is stable when τ<τ00 according to Theorem 1 in [8]. In the presence of diffusion, for example, taking d1=0.008, d2=0.02, which implies (d1+c)(d2+b)=19×10-4>K, we know that the equilibrium E* of system (4) is stable for τ≥0 by Theorem 2.


## 4. The Direction and Stability of Hopf Bifurcation

In this section, we assume the hypotheses of Theorem7 hold and τm=τp=τ. For the case of τm≠τp, which is not our concern in this paper, the calculation of the normal form should follow the method developed in [17]. By using the normal form method in [12] for partial differential equations with time delay, we will investigate the stability of these Hopf bifurcations. For standard notations and classical results on partial functional differential equations, please refer to [12, 13, 17]. More details on techniques for computing the normal form can also be found in recent work [18].Now, normalizing by the time-scalingt→t/τ, then (12)-(13) can be rewritten as
(42)dU(t)dt=τdΔU(t)+L(τ)(Ut)+f(Ut,τ),
where
(43)L(τ)(φ)=τ(a11φ1(0)+a12φ2(-1)a21φ1(-1)+a22φ2(0)),f(φ,τ)=τ(∑i+j≥2‍1i!j!fij(1)φ1i(0)φ2j(-1)∑i+j≥2‍1i!j!fij(2)φ1i(-1)φ2j(0)).In the following, we denote any one of these critical values byτ* at which the characteristic equation (15) has a pair of simply purely imaginary roots ±iω*. Let τ=τ*+α,α∈ℝ, and consider only the case Λ0={-iτ*ω*,iτ*ω*} is the set of eigenvalues on the imaginary axis of the infinitesimal generator associated with the flow of
(44)dU(t)dt=τ*dΔU(t)+L(τ*)(Ut).
Equation (42) is now written as
(45)dU(t)dt=τdΔU(t)+L(τ)(Ut)+F(Ut,α),
where F(φ,α)=αdΔφ(0)+L(α)(φ)+f(φ,τ*+α), for φ∈𝒞. The eigenvalues of τ*dΔ on X are μki=-diτ*k2, i=1,2, k∈ℕ0, with corresponding normalized eigenfunctions βki, where
(46)βk1(x)=(γk(x)0),βk2(x)=(0γk(x)),γk(x)=cos⁡(kx)∥cos⁡(kx)∥2,2,k∈ℕ0.
Let ℬk=span⁡{〈v(·),βki〉βki∣v∈𝒞,i=1,2}. Assume that zt(θ)∈C=C([-1,0],ℝ2) and
(47)ztT(θ)(βk1βk2)∈ℬk.
Then linear PFDE (44) restricted to ℬk is equivalent to the FDE on C([-1,0],ℝ2)(48)z˙(t)=(μk100μk2)z(t)+L(τ*)(zt)
with the characteristic equation given by (15).Suppose that there exists ak∈ℕ0 such that when τ=τ*, (15) for fixed k has a pair of purely imaginary roots ±iω* and all other roots of (15) have negative real parts. Define η(θ)∈BV([-1,0];R) such that
(49)μkψ(0)+L(τ*)ψ=∫-10‍dη(θ)ψ(θ),ψ∈C,
and the adjoint bilinear form on C*×C, C*=C([0,1],ℝ2*),
(50)(ψ(s),ϕ(θ))=ψ(0)ϕ(0)-∫-10‍∫0θ‍ψ(ξ-θ)dη(θ)ϕ(ξ)dξ,000000000000000000000000000000000forψ∈C*,ϕ∈C.
For (48) with fixed k, choose a basis Ψk for the adjoint space P* and a basis Φk for its the eigenspace P as follows:
(51)Φk=(peiω*τ*θ,p-e-iω*τ*θ),Ψk=col(qTe-iω*τ*s,q-Teiω*τ*s)
such that (Φk,Φk)=I2, where I2 is a 2×2 identity matrix. Then we can easily have
(52)p=(p1p2)=(1iω*+d1k2-a11a12eiω*τ*),(53)q=(q1q2)=q1(1iω*+d1k2-a11a21eiω*τ*),
with
(54)q1=(+(iω*+d1k2-a11)2eiω*τ*a12a21)1+2τ*(iω*+d1k2-a11)aa+(iω*+d1k2-a11)2eiω*τ*a12a21)-1.
Following the standard procedure in [12], especially [18], using the decomposition φ(t)=(Φkz)T(βk1βk2)+y, z(t)=(Ψk,(〈φ(·),βk1〉〈φ(·),βk2〉))∈ℝ2, y(t)∈𝒞01∩Ker⁡π=𝒞01∩𝒬:=𝒬1, we decompose (45) as
(55)z˙=Bz+Ψk(0)(〈F((Φkz)T(βk1βk2)+y,α),βk1〉〈F((Φkz)T(βk1βk2)+y,α),βk2〉),ddty=A𝒬1y+(I-π)X0F((Φkz)T(βk1βk2)+y,α),
where here and throughout this section we refer to [12, 18] for results and explanations of several notations involved.Consider the formal Taylor expansion(56)F(v,α)=∑j≥2‍1j!Fj(v,α),
where Fj is the jth Fréchet derivative of F. Then (55) can be written as
(57)z˙=Bz+∑j≥2‍1j!fj1(z,y,α),ddty=A𝒬1y+∑j≥2‍1j!fj2(z,y,α),
where fj1(z,y,α) and fj2(z,y,α) are given by (4.8) in [18]. Then (45) has a normal form on the center manifold of the origin at α=0, written as
(58)z˙=Bz+12g21(z,0,α)+13!g31(z,0,0)+O(α2|z|+α|z|2),
where B=diag⁡{iω*τ*,-iω*τ*} and gj1, j=2,3, are given by (4.9) in [18]. The normal form procedure will show that these terms have the form
(59)12g21(z,0,α)=(Ak1z1αA-k1z2α),
where
(60)Ak1=-k2(d1q1p1+d2q2p2)+iω*qTp;A-k1 is the conjugate of Ak1. Consider the following:
(61)13!g31(z,0,0)=(Ak2z12z2A-k2z1z22),
where
(62)Ak2=i2ω*τ*(ak20ak11-2|ak11|2-13|ak02|2)+12(ak21+12bk21)
with
(63)ak20=τ*∫0π‍γk3(x)dx(b1q1)={τ*π(b1q1),k=0,0,k≠0,b1=f02(1)p22e-2iω*τ*,ak11=τ*∫0π‍γk3(x)dx(b3q1)={τ*π(b3q1),k=0,0,k≠0,b3=f02(1)|p2|2,(64)ak02=τ*∫0π‍γk3(x)dx(b-1q1)={τ*π(b-1q1),k=0,0,k≠0,(65)ak21=τ*∫0π‍γk4(x)dx(q1f03(1)p2|p2|2e-iω*τ*)={τ*π(q1f03(1)p2|p2|2e-iω*τ*),k=0,3τ*2π(q1f03(1)p2|p2|2e-iω*τ*),k≠0,bk21={M0,k=0,M0+22M2k,k≠0,
where, for j=0,2k, Mj=(2τ*/π)q1(f02(1)P2eiω*τ*hj11(2)(-1)+f02(1)P-2e-iω*τ*hj11(2)(-1)), while hk20(θ)=h-k02(θ) and hk20(θ), hk11(θ) are determined by the following equations:
(66)h˙k20(θ)-2iτ*ω*hk20(θ)=Φk(ak20a-k02),h˙k20(0)-L(τ*)(hk20)=τ*ckj(b10),h˙k11(θ)=Φk(2ak112a-k11),h˙k11(0)-L(τ*)(hk11)=τ*ckj(b30),
where
(67)ckj=∫0π‍γk2(x)γj(x)dx={1π,j=k=0,1π,j=0,k≠0,12π,j=2k≠0,0.otherwise.
So the normal form (45) on the center manifold has the form
(68)z˙=Bz+(Ak1z1μA-k1z2μ)+(Ak2z12z2A-k2z1z22)+O(|z|α2+|z4|).
Next we will derive the normal form in the real coordinates. To this end, let z1=w1-iw2, z2=w1+iw2, and then the polar coordinates w1=ρcos⁡ξ, w2=ρsinξ. We finally reach
(69)ρ˙=ιk1αρ+ιk2ρ3+O(α2ρ+|(ρ,α)|4),ξ˙=-ω*τ*+O(|(ρ,α)|);
here ιk1=ReAk1, ιk2=ReAk2. Then, from [19], we know that the number ιk2 tells the bifurcation direction and the stability of bifurcating periodic solution. (i)
Whenιk2<0, it is a supercritical bifurcation and the bifurcating periodic solution is stable.(ii)
Whenιk2>0, it is a subcritical bifurcation and the bifurcating periodic solution is unstable.


## 5. Numerical Simulation and Discussion

In this section, we present some numerical simulations to system (4). These simulations are used to support our theoretical results. Similar to [6, 7, 20], the values of parameters are taken from the published experimental and theoretical results in our simulations where, n∈[2,10], P0∈[40,100], and μ∈[0.01,1]. α=33 and a=4.5 are taken as in [10] for original model (1).Takingd1=0.002, d2=0.02, n=2, b=c=μ=0.03,  and p0=0.4(P0=59.4), then the positive equilibrium E*(m*,p*)=(0.1684,5.6134). From (30) and (32), we obtain the critical value for time delay, τ00≐26.3983. In this case, the parameters satisfy (H-). By Theorem 7, the positive equilibrium E*(0.1684,5.6134) is asymptotically stable for τ=25<τ00. Figure 1 is the numerical simulation of system (4) for τ=25.Figure 1
Numerical simulations of system (4) with d1=0.002, d2=0.02, b=c=0.03, n=2, p0=0.4, and τ=25<τ00. The initial values are m0(x)=0.2+0.1cos⁡x; p0(x)=5.6-0.01cos⁡x. The positive equilibrium E*(0.1684,5.6134) of system (4) is asymptotically stable for ∈[0,τ00).When the delay increasingly crosses through the critical valueτ00≐26.3983, the positive equilibrium E* loses its stability and the Hopf bifurcation occurs. Taking τ=30>τ00, Figure 2 is the numerical simulation results of system (4). It is consistent with the theoretical results.Figure 2
Numerical simulations system (4) with d1=0.002, d2=0.02, b=c=0.03, n=2, p0=0.4, and τ=30>τ00. The initial values are m0(x)=0.2+0.1cos⁡x; p0(x)=5.6-0.01cos⁡x. The positive equilibrium E*(0.1684,5.6134) of system (4) becomes unstable and there exist stable spatially homogeneous periodic solutions.In Figure2, we fix x=1.5708 and the other parameter values are the same as Figure 2. Then we get Figure 3 which shows that the oscillation will sustain when the time delay τ=30 is much greater than the critical value τ00≐26.3983.Figure 3
Sustained oscillation whenτ=30 is much larger than its critical value τ00≐26.3983 with d1=0.002, d2=0.02, b=c=0.03, n=2,p0=0.4 and the initial values are m0(x)=0.2+0.1cos⁡x; and p0(x)=5.6-0.01cos⁡x.Remark 9.
Comparing Figure3 with Figure 3 in [8], we know that the oscillation is still sustained when τ is much greater than its critical value τ00 in the presence of diffusion.



---

*Source: 738682-2014-04-14.xml*