# Stability and Hopf Bifurcation Analysis of an Epidemic Model with Time Delay

Authors: Yue Zhang; Xue Li; Xianghua Zhang; Guisheng Yin
Journal: Computational and Mathematical Methods in Medicine
 (2021)

---

## Abstract

Epidemic models are normally used to describe the spread of infectious diseases. In this paper, we will discuss an epidemic model with time delay. Firstly, the existence of the positive fixed point is proven; and then, the stability and Hopf bifurcation are investigated by analyzing the distribution of the roots of the associated characteristic equations. Thirdly, the theory of normal form and manifold is used to drive an explicit algorithm for determining the direction of Hopf bifurcation and the stability of the bifurcation periodic solutions. Finally, some simulation results are carried out to validate our theoretic analysis.


---

## Body

## 1. Introduction

Today, the serious epidemics, such as SARS and H1N1, are still threatening the life of people continually. Plenty of mathematical models have been proposed to analyze the spread and the control of these diseases [1–7].However, many infectious diseases, for instance, gonorrhea and syphilis, occur and spread amongst the mature, while some epidemics, for example, chickenpox and FMD, only result in infection and death in immature. For this reason, stage structure should be taken into consideration in models. Aliello and Freedman [8] proposed a stage-structured model described by(1)ẋt=αyt−γxt−αe−γτyt−τẏt=αe−γτyt−τ−βy2t,where xt is the immature population density and yt represents the density of the mature population. α, γ, τ, and β are all positive constants. α is the birth rate, and γ is the natural death rate; τ is the time from birth to maturity; β is the death rate of the mature because of the competition with each other.And then, many infectious diseases with sage structure have been built and investigated [9–14]. Xiao and Chen [15] improved (1) by separating the population into mature and immature and supposing that only the immature were susceptible to the infection.Based on the model in [15], supposing that only the mature were susceptible, Jia and Li [16] built a new one as follows:(2)ẋt=αe−γτxt−τ−γxt−βx2t−mxtytẏt=mxtyt−γyt−cyt−gytżt=αxt−γzt−αe−γτxt−τṘt=gyt−γRt,wherext, yt, and Rt are the susceptible, infectious, and recovered mature population densities, respectively; zt denotes the immature population density. All the parameters are positive constants. α, β, γ, and τ are the same as those in (1); m is the transmission coefficient describing the infection between the susceptible and the infectious; c is the death rate because of the epidemic; g is the recovery rate; αe−γτxt−τ denotes the population who were born at t−τ and survive at t.In systems (1) and (2), the time delay was also taken into consideration. Indeed, time delay plays an important role in the epidemic system, making the models more accurate. In recent years, delays have been introduced in more and more epidemic and predator-prey systems [17–19].In this paper, on the basis of (2), we further assume that(1)
Both the susceptible and the infectious have fertility, while in (2), only the susceptible is fertile(2)
For the infectious, there is competition with all the susceptible and the infectious, while for the susceptible, there is only competition between generationsMeanwhile, all the death of the susceptible, the same as that in (1), is only due to the competition. To simplify model (2), we denote γ+c+g=d, and let bxt+yt present the transmission from immature to mature.As a consequence, the new epidemic model could be described as follows:(3)ẋt=bxt+yt−wxt−τ+yt−τxt−mxtytẏt=mxtyt−wxt+ytyt−dytżt=αxt+yt−γzt−bxt+ytṘt=gyt−γRt,where w is the death rate of the mature because of the competition.We can notice thatzt depends on xt and yt and Rt depends on yt; however, xt and yt have nothing to do with zt and Rt. According to Qu and Wei [20], we will mainly focus on xt and yt, that is,(4)ẋt=bxt+yt−wxt−τ+yt−τxt−mxtytẏt=mxtyt−wxt+ytyt−dyt.The rest of the paper is organized as follows. In Section2, we calculate the steady states of system (4) and prove the existence and uniqueness of the positive equilibrium in particular. And then, the stability of the two nonzero equilibria and the existence of the Hopf bifurcation are investigated in Sections 3 and 4, respectively. In Section 5, the direction and stability of the Hopf bifurcation at the positive equilibrium are studied by using the center manifold theorem and the normal form theory [21]. And in the last section, some numerical simulations are carried out to validate the theoretical analysis.


## 2. The Existence and Uniqueness of the Positive Equilibrium of the Model

In this section, we discuss the existence of the equilibria of (4) and the positive one in particular.The equilibria are the solutions of the equations (5),(5)bx+y−wx+yx−mxy=0mx−wx+y−dy=0.Clearly,E10,0 and E2b/w,0 are two equilibria of (4).In the following, we will focus on the existence of the positive equilibrium.Theorem 1.
Ifbm−w>dw, (4) has one positive equilibrium E3x∗,y∗, where(6)x∗=B+∆2m2,y∗=m−wB+∆/2m2w−d/w.Proof.
Positive equilibrium is the positive solution of the equations (7),(7)bx+y−wx+yx−mxy=0mx−wx+y−d=0.
From the second equation of (7), we have(8)wy=m−wx−d.
Taking (8) into the first equation of (7), we can obtain(9)m2x2−bm+dm+dwx+bd=0,
which leads to(10)∆=B2−4m2bd=bm+dw−dm2+4d2mw>0,
where(11)B=bm+dm+dw.
Together with(12)x1+x2=bm+dm+dw/m2>0,andx1x2=bd/m2>0,
we can know that both of the two solutions of (9) are positive,
where(13)x1=B−∆2m2>0andx2=B+∆2m2>0.
Ifx=B−∆/2m2,
then,(14)wy=m−wx−d=m−wB−∆2m2−d<m−wB−bm+dw−dm2m2−d≤m−wB−bm+dw−dm2m2−d=−dwm<0.
So,x=B−∆/2m2 is dropped.
Ifx=B+∆/2m2,
then,(15)wy=m−wB+∆2m2−d=m−w4m2bd2m2B−∆−d=d2bm−wB−∆−1=dB−∆2bm−w−B−∆.
w>0, d>0, and B−∆=∆+4m2bd−∆>0, so y>0⟺2bm−w−B−∆>0,(16)⟺∆>B2−4bm−wB+4b2m−w2⟺m2d−m−wbm−bw−bm−dm−dw<0⟺bm−w>dw∴bm−w>dw⟺y>0.
Then, takingx=B+∆/2m2≜x∗ into (8), we have(17)y=m−wB+∆2m2w−dw≜y∗.
Therefore, ifbm−w>dw, (4) has the unique positive equilibrium E3x∗,y∗.☐


## 3. Stability Analysis of the EquilibriumE2b/w,0

In this section, we analyze the stability of the equilibriumE2b/w,0.For convenience, the new variablesut=xt−b/w and vt=yt are introduced, and then, around E2b/w,0, the system (4) could be linearized as (18):(18)u̇t=−but−τ−bvt−τ+bvt−bmwvtv̇t=bmwvt−dvt−bvt,whose characteristic equation is given by(19)λ+be−λτλ+b+d−bmw=0,from which, we can get that(20)λ=−b+d−bmw,or(21)λ+be−λτ=0.Obviously, ifbm−w>dw,then(22)λ=−b+d−bmw>0,which implies that the equilibriumE2 is unstable.Ifbm−w<dw,then(23)λ=−b+d−bm/w<0.As a consequence, we will discuss the other roots of (19), that is, the roots of (21), under the condition bm−w<dw.Forτ=0, equation (21) becomes(24)λ+b=0,whose root is(25)λ=−b<0.Forτ>0, if iωω>0 is a root of (21), then(26)iω+bcosωτ−ibsinωτ=0.(27) can be obtained by separating the real and the imaginary parts,(27)bcosωτ=0bsinωτ=ω,which leads to(28)ω2=b2,from which, we can get the unique positive root(29)ω0=b.Let(30)τj=π4j+12ω0j=0,1,2,⋯.Then, whenτ=τj, (21) has a pair of purely imaginary roots ±iω0.Suppose(31)λτ=ατ+iωτ,which is the root of (21) such that(32)ατj=0,andωτj=ω0.To investigate the distribution ofλτ, we will discuss the trend of λτ at τ=τj.Substitutingλτ into (21) and taking the derivative with respect to τ, we can get(33)dλdτ−be−λτλ+τdλdτ=0,which yields,(34)dλdτ−1=eλτbλ−τλ.Together with (27), we have(35)signRedλdττ=τj−1=signReeλτbλτ=τj=signRecosω0τ+isinω0τibω0=signsinω0τbω0=signbsinω0τb2ω0=signω0b2ω0=sign1b2>0,which means that when undergoesτ=τj, λτ will add a pair of roots with positive real parts. That is, with the increase of τ, the number of roots with positive real part is increasing, leading to the change of the stability of the system (4).Therefore, the distribution of the roots of (21) could be obtained.Lemma 2.
Letω0 and τj (j=0,1,2,⋯) be defined by (29) and (30), respectively.(1)
Ifbm−w>dw, then (19) has at least one positive root(2)
Ifbm−w<dw, and τ=0, then both roots of (19) are negative(3)
Ifbm−w<dw, and τ>0, then (19) has a pair of simple imaginary roots ±iω0 at τ=τj; furthermore, if τ<τ0, then all the roots of (19) have negative real parts; if τ∈τj,τj+1, (19) has 2j+1 roots with positive real partsTogether with condition (35), the Hopf bifurcation theorem [21], and Lemma 2, the following theorem could be stated.Theorem 3.
Letτj (j=0,1,2,⋯) be defined by (30), then we have(1)
Ifbm−w>dw, then the equilibrium E2b/w,0 of (4) is unstable(2)
Ifbm−w<dw, then the equilibrium E2b/w,0 is asymptotically stable when τ∈0,τ0, and it is unstable when τ>τ0(3)
Ifbm−w<dw, then system (4) undergoes a Hopf bifurcation at the equilibrium E2b/w,0 for τ=τj


## 4. Stability Analysis of Positive EquilibriumE3x∗,y∗

In this section, we analyze the stability of the positive equilibriumE3x∗,y∗.For convenience, the new variablesut=xt−x∗ and vt=yt−y∗ are introduced, and then, around E3x∗,y∗, the system (4) could be linearized as (36):(36)u̇t=b+d−mx∗−my∗ut−wx∗ut−τ+b−mx∗vt−wx∗vt−τv̇t=m−wy∗ut−wy∗vt,whose characteristic equation is given by(37)λ2+wy∗−kλ+λ+my∗se−λτ−kwy∗−m−wy∗n=0,where(38)k=b+d−mx∗−my∗,n=b−mx∗,s=wx∗.Forτ=0, equation (37) becomes(39)λ2+wy∗−k+sλ+msy∗−kwy∗−m−wy∗n=0.Firstly, computingλ1λ2, we have(40)λ1λ2=msy∗−kwy∗−m−wy∗n=y∗mwx∗−b+dw+mx∗+y∗w−m−wb+m−wmx∗=y∗m2x∗+mwx∗+mmx∗−wx∗−d−dw+bm=y∗Δ,whereΔ is the same as that in (10).So,(41)λ1λ2=msy∗−kwy∗−m−wy∗n>0,which implies that the real parts ofλ1 and λ2 have the same signs.Then, (λ1+λ2) is calculated:(42)λ1+λ2=−wy∗−k+s=b+d−mx∗+y∗y∗=b+mx∗−wx∗+y∗−mx∗−wx∗+y∗−my∗=−my∗2x∗+y∗−wx∗+y∗<0.Together with (41), we can get that both the real parts of the two roots of (39) are negative.Forτ>0, equation (37) can be rewritten as(43)λ2+rλ+sλ+ce−λτ+p=0,where(44)r=wy∗−k,s=wx∗,c=msy∗,p=−kwy∗−m−wy∗n.Ifiωω>0 is a root of (37), then(45)−ω2+irω+isωcosωτ+swsinωτ+ccosωτ−icsinωτ+p=0.Separating the real and the imaginary parts, we have(46)ccosωτ+sωsinωτ=ω2−pcsinωτ−sωcosωτ=rω,which leads to(47)ω4−2p−r2+s2ω2+p2−c2=0.Letz=ω2>0, and then, (47) can be rewritten as(48)z2−2p−r2+s2z+p2−c2=0.Firstly, computingz1z2, we get(49)z1z2=p2−c2=p+cp−c,where(50)p−c=msy∗−kwy∗−m−wy∗n>0,has been proved in (41).Then, we will calculate(51)p+c=−msy∗−kwy∗−m−wy∗n=y∗−mwx∗−b+dw+mx∗+y∗w−m−wb+m−wmx∗=y∗m2x∗−mwx∗+mwy∗−dw+bm=y∗m−wmB+Δ−B=y∗mm−wΔ−wB.If H(4-1):m−wΔ−wB<0,then(52)p+c<0,and then,(53)z1z2<0,which implies that (48) has one unique positive solution z0=ω02,where(54)ω0=2p−r2+s2+2p−r2+s22−4p2−c221/2.If H(4-2):m−wΔ−wB>0,then(55)p+c>0,and then,(56)z1z2>0,Let(57)∆1=2p−r2+s22−4p2−c2.If∆1<0, then (48) has no real roots.If∆1>0 and z1+z2=2p−r2+s2<0, then (48) has two negative roots, and there no positive ω for (47);If∆1>0 and z1+z2=2p−r2+s2>0,then (48) has two positive roots, and there are two positive ω for (47), which are(58)ω±=2p−r2+s2±2p−r2+s22−4p2−c221/2.Lemma 4.
(1)
Ifm−wΔ−wB<0, then (47) has one positive root ω0(2)
Whenm−wΔ−wB>0, if ∆1<0, or z1+z2=2p−r2+s2<0, then (47) has no positive roots(3)
Whenm−wΔ−wB>0, if ∆1>0 and z1+z2=2p−r2+s2>0, then (47) has two positive rootsIn the following, we will discuss the expression ofτj.Lemma 5.
If(59)cosat=fasinat=ga,
wherea>0, t>0, then(1)
Iffa>0,ga>0, then(60)at=arccosfa+2jπ,orat=arcsinga+2jπ(2)
Iffa<0,ga>0, then(61)at=arccosfa+2jπ,orat=π−arcsinga+2jπ(3)
Iffa>0,ga<0, then(62)at=2π−arccosfa+2jπ,orat=2π+arcsinga+2jπ(4)
Iffa<0,ga<0, then(63)at=2π−arccosfa+2jπ,orat=π−arcsinga+2jπ
In conclusion,(1)
Ifga>0, then at=arccosfa+2jπ(2)
Ifga<0, then at=2π−arccosfa+2jπ,where j=0,1,2,⋯According to (46), we have(64)cosωτ=cω2−p−srω2c2+s2ω2=fωsinωτ=sωω2−p+crωc2+s2ω2=gω.Ifm−wΔ−wB<0, substituting ω0 defined in (54) into (64), we can get fω0 and gω0. Together with Lemma 5, the expression of τj could be obtained.If(65)gω0>0,then(66)τj=1ω0arccoscω02−p−srω02c2+s2ω02+2πj,j=0,1,2,⋯.If(67)gω0<0,then(68)τj=1ω02π−arccoscω02−p−srω02c2+s2ω02+2πj,j=0,1,2,⋯.That is, whenτ=τj, the characteristic equation (37) has a pair of purely imaginary roots ±iω0.Supposeλτ=ατ+iωτ is the root of (37), and then, we have(69)ατj=0,andωτj=ω0.To investigate the distribution of theλτ, we will discuss the trend of λτ at τ=τj.Substitutingλτ into (37) and taking the derivative with respect to τ, we can get(70)2λ+r+se−λτ−e−λτsλ+cτdλdτ=λe−λτsλ+c,which leads to(71)dλdτ−1=eλτsλ+c+sλsλ+c−τλ.Together with (46) and (54), we have(72)signRedλdττ=τj−1=signReeλτsλ+c+sλsλ+cτ=τj=signRecosω0τ+isinω0τ2iω0+r+siω0isω0+c=signc2ω0cosω0τ+rsinω0τ−sω0rcosω0τ−2ω0sinω0τ+s=sign2ω0ω02−p+rrω0=sign2p−r2+s2+2p−r2+s22−4p2−c2−2p−r2+s2=sign2p−r2+s22−4p2−c2>0,which means that when undergoesτ=τj, λτ will add a pair of roots with positive real parts. That is, with the increase of τ, the number of roots with positive real part is increasing, leading to the change of the stability of the system (4).Ifm−wΔ−wB>0, ∆1>0 and z1+z2=2p−r2+s2>0, then fω0 and gω0 could be calculated by substituting ω± defined in (58) into (64). According to Lemma 5, we can get the expression of τj±:If(73)gω±>0,then(74)τj±=1ω±arccoscω±2−p−srω±2c2+s2ω±2+2πj,j=0,1,2,⋯.If(75)gω±<0,then(76)τj±=1ω±2πj+1−arccoscω±2−p−srω±2c2+s2ω±2,j=0,1,2,⋯.That is, whenτ=τj±, the characteristic equation (37) has a pair of purely imaginary roots ±iω±.Letλτ=ατ+iωτ be the root of (37), satisfying(77)ατj±=0,andωτj±=ω±.To investigate the distribution of theλτ, we will discuss the trend of λτ at τ=τj±.Using the same method, we have(78)signRedλdττ=τj±−1=sign2ω±2−2p−r2+s2=sign2p−r2+s2±2p−r2+s22−4p2−c2−2p−r2+s2=sign±2p−r2+s22−4p2−c2=sign±∆1.This implies that(79)dλdττ=τj+>0,and(80)dλdττ=τj−<0,which means that when undergoesτ=τj+, λτ will add a pair of roots with positive real parts, while undergoes τ=τj−, λτ will lose a pair of roots with positive real parts; if τj−>τj+1+, then the characteristic equation (37) must have roots with positive real parts for τ>τj+1+.In conclusion, the distribution of the roots of (37) could be obtained.Lemma 6.
Letω0 be defined by (54), and τj (j=0,1,2,⋯) be defined by (66) or (68), and ω± be defined by (58), and τj± (j=0,1,2,⋯) be defined by (74) or (76), respectively.(1)
Whenm−wΔ−wB>0, if ∆1<0, or z1+z2=2p−r2+s2<0, then all the roots of (37) are with negative real parts(2)
Whenm−wΔ−wB<0, (37) has a pair of simple imaginary roots ±iω0 at τ=τj; furthermore, if τ∈0,τ0, then all the roots of (37) are with negative real parts; if τ∈τj,τj+1, then (37) has 2j+1 roots with positive real parts(3)
Whenm−wΔ−wB>0, if ∆1>0, and z1+z2=2p−r2+s2>0, then (37) has a pair of simple imaginary roots ±iω± at τ=τj±; if τ∈0,τ0+ or τ∈τj−,τj+1+, then all the roots of (37) are with negative real parts; if τ∈τj+1+,τj+1−, then (37) has a pair of roots with positive real parts; if τh−>τh+1+, for τ>τh+1+, (37) must have roots with positive real partsTogether with conditions (72) and (78), the Hopf bifurcation theorem [21], and Lemma 6, the following theorem could be stated.Theorem 7.
Letτj (j=0,1,2,⋯) be defined by (66) or (68), and τj± (j=0,1,2,⋯) be defined by (74) or (76), respectively, and then we have(1)
Whenm−wΔ−wB>0, if ∆1<0, or z1+z2=2p−r2+s2<0, the equilibrium E3x∗,y∗ of (4) is asymptotically stable(2)
Whenm−wΔ−wB<0, then the equilibrium E3x∗,y∗ is asymptotically stable when τ∈0,τ0, and it is unstable when τ>τ0. System (4) undergoes a Hopf bifurcation at the equilibrium E3x∗,y∗ for τ=τj(3)
Whenm−wΔ−wB>0, if ∆1>0, and z1+z2=2p−r2+s2>0, then there is a positive h, such that when τ∈0,τ0+∪τ0−,τ1+∪τ1−,τ2+∪⋯∪τh−1−,τh+, the equilibrium E3x∗,y∗ is asymptotically stable, and it is unstable when τ∈τ0+,τ0−∪τ1+,τ1−∪τ2+,τ2−∪⋯∪τh−1+,τh−1−∪τh+,+∞. We call that system (4) undergoes h+1 switchesSystem (4) undergoes a Hopf bifurcation at the equilibrium E3x∗,y∗ for τ=τj±.


## 5. The Direction and Stability of Hopf Bifurcation atE3x∗,y∗

In the previous section, we have already gotten some conditions making that the system (4) undergoes a Hopf bifurcation at the positive equilibrium E3x∗,y∗ when τ=τj±, j=0,1,2,⋯. In this section, under the conditions in Theorem 7, the direction of Hopf bifurcation and stability of the periodic solutions from E3 will be investigated by using the center manifold and normal form theories [21].Without loss of the generality, letτ¯ be the critical value of τ=τj+τj−, j=0,1,2,⋯.For convenience, letτ=τ¯+ρ, ρ∈R, ut=xτt−x∗, vt=yτt−y∗, and then system (4) undergoes Hopf bifurcation at ρ=0; with τ normalized by the time scaling t⟶t/τ, (4) could be rewritten as(81)u̇t=τbut+vt+x∗+y∗−wut−1+vt−1+x∗+y∗ut+x∗−mut+x∗vt+y∗v̇t=τmut+x∗vt+y∗−dvt+y∗−wut+x∗vt+y∗.Choose the space asC=C−1,0,R2; for any Φ=Φ1,Φ2∈C, let(82)LρΦ=τ¯+ρmΦ10−sΦ1−1+nΦ20−sΦ2−1m−wy∗Φ10−wy∗Φ20,and(83)fτ¯,ρ=τ¯+ρ−wΦ10Φ1−1+Φ2−1−mΦ10Φ20m−wΦ10Φ20−wΦ220,wheres andn are the same as those in (37).According to Riesz’s representation theorem, there is a2×2 matrix ηθ,ρ (θ∈−1,0), whose components are bounded variation functions such that(84)LρΦ=∫−10dηθ,ρΦθforθ∈C.In fact,ηθ,ρ could be chosen as(85)ηθ,ρ=τ¯+ρmnm−wwy∗δθ+τ¯+ρ−s−s00δθ+1,where(86)δθ=0,θ=01,θ≠0.For anyΦ∈C1−1,0,R2, we define(87)AρΦθ=dΦθdθ,θ∈−1,0∫−10dησ,ρΦσ,θ=0,and(88)RρΦ=0,θ∈−1,0fρ,Φm,θ=0.And then, system (4) could be translated into(89)ẋt=Aρxt+Rρxt,where(90)xt=xt+ρ=ut+ρvt+ρ.ForΨ∈C1−1,0,R2∗, define(91)A∗Ψs=−dΨsds,s∈−1,0∫−10dηt,0Ψ−t,s=0,and a bilinear form(92)Ψs,Φθ=Ψ¯0Φ0−∫−10∫ξ=0θΨ¯ξ−θdηθΨξdξ,where(93)ηθ=ηθ,0.According to (87) and (91), we can get that A0 and A∗ are adjoint operators.From the analysis in Section4, we know that ±iωτ¯ are a pair of eigenvalues of A0 and also eigenvalues of A∗, where ω is ω+ or ω− defined in (58).It is easy to verify that vectorsqθ and q∗s are the eigenvalues of A0 and A∗ corresponding to the eigenvalues iωτ¯ and −iωτ¯, where(94)qθ=iω+wy∗m−wy∗,1Teiωτ¯θ,q∗s=B¯1,−iω−m+seiωτ¯m−wy∗eiωτ¯s.For convenience, let(95)M¯=−iω−m+seiωτ¯m−wy∗,N=iω+wy∗m−wy∗,then(96)q∗s,qθ=B1,MN1−∫−10∫ε=0θ1,Me−iωτ¯ε−θdηθN1eiωτ¯εdε=BM+N−B∫−101,Mθeiωτ¯θdηθN1=BM+N−B1,M∫−10dηθθNeiωτ¯θθeiωτ¯θ=BM+N−B1,Mτ¯sN+1e−iωτ¯0=BM+N−Bτ¯sN+1e−iωτ¯.We choose(97)B=M+N−τ¯sN+1e−iωτ¯−1,and then,(98)q∗s,qθ=1,qθ,q∗s=0.Using the same method as that in [21], the center manifold C0 at ρ=0 is first computed. Suppose that xt is the solution of (81) when ρ=0, and define(99)zt=q∗,xtθ,Wt,θ=xtθ−2Reztqθ.Then, on the center manifoldC0, we can get that(100)Wt,θ=Wzt,z¯t,θ,where(101)Wzt,z¯t,θ=W20θz22+W11θzz¯+W02θz¯22+W30θz32+⋯.In the directionq∗ and q¯∗, z and z¯ are local coordinates for center manifold C0. It is not difficult to note that when xt is real, W is real. And the real solutions are considered only.For any solutionxt∈C0, since ρ=0, we can get that(102)żt=q∗,xṫ=iωτ¯zt+q¯∗0f0,Wzt,z¯t,θ+2Reztq0≝iωτ¯zt+q¯∗0f0z,z¯,where(103)f0z,z¯=fz2z22+fz¯2z¯22+fzz¯zz¯+fz2z¯z2z¯2+⋯.We rewrite(104)żt=iωτ¯zt+gz,z¯,with(105)gz,z¯=q¯∗0f0z,z¯=g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯.Together with (83), we have(106)gz,z¯=q¯∗0f0,Wzt,z¯t,θ=q¯∗0f0,xt=1,MBτ¯−wW10+zN+z¯N¯W1−1+zNe−iωτ¯+z¯N¯eiωτ¯+W2−1+ze−iωτ¯+z¯eiωτ¯−mW10+zN+z¯N¯W20+z+z¯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯m−wW10+zN+z¯N¯W20+z+z¯−wW20+z+z¯2=−2Bτ¯wN2e−iωτ¯+Ne−iωτ¯+mNz22−Bτ¯wNN¯+1eiωτ¯+wN¯N+1e−iωτ¯+mN+N¯zz¯−2Bτ¯wN¯N+1eiωτ¯+mN¯z¯22−Bτ¯wW2010N¯+1eiωτ¯+2W1110N+1e−iωτ¯+2NW112−1+W111−1+mW2010+2W1110+2NW1120+N¯W2020z2z¯2+2Bτ¯Mm−wN−wz22+Bτ¯MN+N¯m−w−2wzz¯+2Bτ¯Mm−wN¯−wz¯22+Bτ¯Mm−wW2010+2W1110+2NW1120−2w2W1120+W2020z2z¯2.Comparing with (105), we get(107)g20=−2Bτ¯wN2e−iωτ¯+Ne−iωτ¯+mN+2Bτ¯Mm−wN−w,g11=−Bτ¯wNN¯+1eiωτ¯+wN¯N+1e−iωτ¯+mN+N¯+Bτ¯MN+N¯m−w−2wg02=−2Bτ¯wN¯N+1eiωτ¯+mN¯+2Bτ¯Mm−wN¯−wg21=−Bτ¯wW2010N¯+1eiωτ¯+2W1110N+1e−iωτ¯+2NW112−1+W111−1+mW2010+2W1110+2NW1120+N¯W2020+Bτ¯Mm−wW2010+2W1110+2NW1120−2w2W1120+W2020.All the numbers in the expressions ofg20, g11, and g02 are known; however, W20 and W11 in g21 are unknown, which will be computed in the following.From (89) and (99),(108)Ẇ=xṫ−żq−z¯̇q¯=AW−2Req¯∗0f0qθ,θ∈−1,0AW−2Req¯∗0f0q0+f0,θ=0≝AW+Hzt,z¯t,θ,where(109)Hzt,z¯t,θ=H20θz22+H11θzz¯+H02θz¯22+⋯.Whenθ∈−1,0,(110)Hzt,z¯t,θ=−q¯∗0f0qθ−q∗0f0¯q¯θ=−gz,z¯qθ−g¯z,z¯q¯θ.Comparing the coefficients in the expansions of (109) and (110), we have(111)H20θ=−g20qθ−g¯02q¯θH11θ=−g11qθ−g¯11q¯θ.Taking (5-18) into (5-17), we obtain(112)A−2iωτ¯W20θ=−H20θAW11θ=−H11θ⋯.By (111) and (112), the following can be gotten(113)Ẇ20θ=2iωτ¯W20θ+g20qθ+g¯20q¯θ,Ẇ11θ=g11qθ+g¯11q¯θ.Because(114)qθ=q0eiωτ¯θ,by integration, we have(115)W20θ=ig20ωτ¯q0eiωτ¯θ+ig¯203ωτ¯q¯0e−iωτ¯θ+E1e2iωτ¯θ,and(116)W11θ=−ig11ωτ¯q0eiωτ¯θ+ig¯11ωτ¯q¯0e−iωτ¯θ+E2,whereE1 and E2 are unknown.From (112) and the definition of A, we obtain(117)∫−10dηθW20θ=2iωτ¯W20θ−H20θ,and(118)∫−10dηθW11θ=−H11θ.Then, together with (108), we get(119)H200=−g20q0−g¯02q¯0+2−wNNe−iωτ¯+1−mNm−wN−wH110=−g11q0−g¯11q¯0+2−wReNNe−iωτ¯+1−mReNm−wReN−w.Substituting (115) and (119) into (117), we have(120)2iωτ¯I−∫−10eiωτ¯θdηθE1=2−wNNe−iωτ¯+1−mNm−wN−w.That is(121)2iω−m+se−2iωτ¯−n+se−2iωτ¯−m−wy∗2iω+wy∗E1=2−wNNe−iωτ¯+1−mNm−wN−w,from which,E1 can be determined.By (115), (117), and (119), we have(122)∫−10dηθE2=2−wReNNe−iωτ¯+1−mReNm−wReN−w.That is(123)−m+s−n+s−m−wy∗wy∗E2=2−wReNNe−iωτ¯+1−mReNm−wReN−w,from which,E2 can be determined.SubstitutingE1 and E2 into (115) and (116), W20 and W11 could be obtained; furthermore, g21 can be calculated.Then, the important parameters can be obtained:(124)C10=i2ωτ¯g20g11−2g112−13g022+g212,μ2=−ReC10Reλ′τ¯,β2=2ReC10,T2=−ImC10+μ2Imλ′τ¯ωτ¯,which determine the quantities of the bifurcation of system (4) at τ=τ¯, where μ2 determines the direction of the bifurcation: Hopf bifurcation is subcritical (supercritical) if μ2<0 (μ2>0) and the bifurcating periodic solutions exist for τ<τ¯ (τ>τ¯); β2 determines the stability of the bifurcating periodic solutions, which are stable (unstable) if β2<0 (β2>0); T2 determines the period of the periodic solutions, which decreases (increases) if T2<0 (T2>0).From (78) in Section 4, we know that Reλ′τ¯>0 when τ=τj+, and Reλ′τ¯<0 when τ=τj−.Theorem 8.
Ifm−wΔ−wB<0, then the Hopf bifurcation of system (4) at positive equilibrium E3x∗,y∗ when τ=τj is supercritical (subcritical) and the bifurcating periodic solutions on the manifold are stable (unstable) if ReC10<0 (ReC10>0).Theorem 9.
Ifm−wΔ−wB>0, ∆1>0, and z1+z2=2p−r2+s2>0, then when τ=τj+, the Hopf bifurcation of system (4) at the positive equilibrium E3x∗,y∗ is supercritical (subcritical) and the bifurcating periodic solutions on the manifold are stable (unstable) if ReC10<0 (ReC10>0);
whenτ=τj−, the Hopf bifurcation of system (4) at the positive equilibrium E3x∗,y∗ is subcritical (supercritical) and the bifurcating periodic solutions on the manifold are stable (unstable) if ReC10<0 (ReC10>0).


## 6. Numerical Simulations

In this section, some numerical simulations are carried out to support our theoretical analysis.There are so many different cases that only the most particular one in (3) of Theorem7 is considered in this section. The coefficients are chosen as follows: b=0.4, d=0.8, w=1, andm=6. Then, the conditions of (3) in Theorem 7 are satisfied, wherem−wΔ−wB≈13.166>0,∆1≈0.0110067>0,and2p−r2+s2≈0.319612>0.By direct calculation, we have(125)x∗≈0.169906,y∗≈0.049528;and(126)τ0+≈4.25918,τ0−≈9.66231,τ1+≈17.8969,τ1−≈28.8393,τ2+≈31.5347,τ2−≈48.0163,τ3+≈45.1725⋯,whereτ2−>τ3+.From Theorem7, we know that the positive equilibrium E30.169906,0.049528 should be asymptotically stable when τ∈0,τ0+∪τ0−,τ1+∪τ1−,τ2+, and it is unstable when τ∈τ0+,τ0−∪τ1+,τ1−∪τ2+,+∞. The system (4) undergoes 3 switches.All the simulation results forτ in the six different intervals are in consistent with the theoretical analysis, which are shown in Figures 1–6.Figure 1
The dynamic behaviors of system (4) with τ=3.8∈0,τ0+. (a) and (b) are the waveform plot and phase for system (4), respectively. E3x∗,y∗ is asymptotically stable and the initial condition is 0.16,0.1.
(a)(b)Figure 2
The dynamic behaviors of system (4) with τ=4.4∈τ0+,τ0−. (a) and (b) are the waveform plot and phase for system (4), respectively. E3x∗,y∗ is unstable with the emergence of oscillation and the initial condition is 0.18,0.12.
(a)(b)Figure 3
The dynamic behaviors of system (4) with τ=11∈τ0−,τ1+. (a) and (b) are the waveform plot and phase for system (4), respectively. E3x∗,y∗ is asymptotically stable and the initial condition is 0.1,0.06.
(a)(b)Figure 4
The dynamic behaviors of system (4) with τ=23∈τ1+,τ1−. (a) and (b) are the waveform plot and phase for system (4), respectively. E3x∗,y∗ is unstable with the emergence of oscillation and the initial condition is 0.2,0.2.
(a)(b)Figure 5
The dynamic behaviors of system (4) with τ=30∈τ1−,τ2+. (a) and (b) are the waveform plot and phase for system (4), respectively. E3x∗,y∗ is asymptotically stable and the initial condition is 0.15,0.02.
(a)(b)Figure 6
The dynamic behaviors of system (4) with τ=36∈τ2+,+∞. (a) and (b) are the waveform plot and phase for system (4), respectively. E3x∗,y∗ is unstable with the emergence of oscillation and the initial condition is 0.15,0.14.
(a)(b)



---

*Source: 1895764-2021-07-02.xml*