# Bifurcations and Structures of the Parameter Space of a Discrete-Time SIS Epidemic Model

Authors: Yun Liu; Xijuan Liu
Journal: Journal of Mathematics
 (2022)

---

## Abstract

The dynamics of a discrete-time SIS epidemic model are reported in this paper. Three types of codimension one bifurcation, namely, transcritical, flip, Neimark–Sacker (N-S) bifurcations, and their intersection codimension two bifurcations including 1 : 2, 1 : 3, and 1 : 4 resonances are discussed. The necessary and sufficient conditions for detecting these types of bifurcation are derived using algebraic criterion methods. Numerical simulations are conducted not only to illustrate analytical results but also to exhibit complex behaviors which include period-doubling bifurcation in period−2,−4,−8,−16 orbits, invariant closed cycles, and attracting chaotic sets. Especially, here we investigate the parameter space of the discrete model. We also investigate the organization of typical periodic structures embedded in a quasiperiodic region. We identify period-adding, Farey sequence of periodic structures embedded in this quasiperiodic region.


---

## Body

## 1. Introduction

Throughout the history of the world, epidemic diseases have posed a lasting threat to humans and society. For instance, tuberculosis, Ebola virus, Zika virus, influenza, etc. may bring many disasters to public security and individual health. Therefore, it has become an urgent and realistic problem to study the pathogenesis, transmission mechanism, development trend, and control strategy of infectious diseases. It is of great practical significance to analyze the spread of diseases and propose appropriate control strategies by using mathematical models and their dynamical properties [1–4]. The insights gained from mathematical modeling in epidemiology are important in the fight against epidemics.In general, the dynamics of epidemics are a highly transited topic of investigation from the epidemiological point of view. There are various families of models depending on the assumption of the mechanisms of propagation. For example, there are susceptible-infected-susceptible (SIS) systems that have been used to model nonlinear incidence rates and double epidemic hypotheses [5], susceptible-infected-recovered (SIR) models which are employed to consider the awareness of the presence of a disease [6], susceptible-exposed-infected-recovered (SEIR) systems that predict the propagation of epidemics with nonlocal reaction functions [7], susceptible-exposed-infected-quarantined-recovered (SEIQR) systems used to account for adjusted incidence and imperfect vaccinations [8], and so on. Generally, in the theory of epidemic dynamical models, there are two kinds of mathematical models: the continuous-time models described by differential equations, and the discrete-time models described by difference equations. In recent years, more and more attention has been paid to the discrete-time epidemic models [9–12].As is well known, difference equations have many applications in epidemic models [13–15]. In real life, data collection of disease is often based on a certain interval, and discrete models are more convenient for data simulation and comparison. On the other hand, a difference equation is the discretization of the continuous model, which makes it easier to handle for the numerical approach. Especially, discrete models exhibit more complex dynamical behaviors than the corresponding continuous models. With these advantages in mind, the problem of discrete epidemic models is of great research interest [16–18]. Based on the abovementioned reasons, we decide to consider a discrete-time epidemic model, which will be meaningful work. In this paper, we will study a discrete-time SIS epidemic model given by Hu et al. [19].(1)Sn+1=Sn+hA−λSnInSn+In−d1Sn+rIn,In+1=In+hλSnInSn+In−d2+rIn,where Sn and In denote the numbers of susceptible and infective individuals at time t, respectively; A is the recruitment rate of the population; d1 is the natural death rate of the population; d2 is the death rate of infective individuals which includes the natural death rate and the disease-related death rate; r stands for the recovery rate of the infective individuals; and λ is the standard incidence rate. It is assumed that initial values S0,I0>0 and all the parameters are positive constants. And, h is the step size.For the sake of simplicity, we rewrite system (1) in the following form:(2)xy⟶x+hA−λxyx+y−d1x+ryy+hλxyx+y−d2+ry.The objective is to study systematically the dynamics of model (2). Some of the codimension one bifurcations of the above model are investigated in [19]. Here, all of the bifurcations including codimension one and codimension two along with the normal form coefficients are provided; for more details, see [20–22].The remainders of the paper are organized as follows: In Section2, we present the existence and local stability of fixed points for map (2). In Section 3, we show that there exist some values of parameters such that map (2) possesses codimension one bifurcation. In Section 4, we extract the codimension two bifurcations of the model. In Section 5, we provide some numerical simulations, which not only validate the effectiveness of the theoretical analysis but also illustrate the mode-locking structures in the quasiperiodic rotation motions. Finally, we conclude this paper with comments and discussions.


## 2. Existence and Stability of Fixed Points

In this section, we discuss the existence of the positive fixed points of model (2). The basic reproductive rate of model (2) is R0=λ/d2+r, which is the average number of secondary infections generated by an initial population of infected individuals over their lifetimes. If R0<1, model (2) has a unique fixed point E0A/d1,0. If R0>1, model (2) has two fixed points: E0A/d1,0 and E1x1,y1, where(3)x1=Ad2+rd1r+d2+d2λ−r−d2,y1=Aλ−r−d2d1r+d2+d2λ−r−d2.Let us consider model (2), given by the map(4)xy⟶x+hA−λxyx+y−d1x+ryy+hλxyx+y−d2+ry=Φx,y,ϑ,where ϑ=h,A,d1,d2,r,λ. Its Jacobian matrix J=ΦX, where X=x,y is given by(5)Jx,y=1−hd1+λy2x+y2hr−λx2x+y2hλy2x+y21+hλx2x+y2−r+d2.And, the multilinear functions of (2) used in its Taylor expansion can be expressed as follows:(6)BiΓ,Σ=∑j,k=12∂2Φix,y∂ξj∂ξkrjσk,CiΓ,Σ,ϒ=∑j,k,l=12∂3Φix,y∂ξj∂ξk∂ξLrjσkvl,where(7)ξ=x,y,Γ=r1,r2T,Σ=σ1,σ2T,ϒ=v1,v2T.Thus,(8)BΓ,Σ=h1,h2T,CΓ,Σ,ϒ=g1,g2T,where(9)h1=2hλx+y3y2r1σ1−xyr1σ2+r2σ1+x2r2σ2,h2=−h1,g1=2hλx+y4−3y2r1σ1v1+2xy−y2r1σ1v2+r2σ1v1+r1σ2v1−3x2r2σ2v2+2xy−x2r2σ2v1+r1σ2v2+r2σ1v2,g2=−g1.


## 3. Codimension One Bifurcation Analysis

The following theorems provide the codimension one bifurcation analysis for fixed pointsE0 and E1.


### 3.1. Codimension One Bifurcation atE0

Theorem 1.
(1)
There is a transcritical bifurcation ofE0 at λ=r+d2 when hd1≠2(2)
There is a transcritical bifurcation ofE0 at r=λ−d2 when hd1≠2(3)
There is a transcritical bifurcation ofE0 at d2=λ−r when hd1≠2Proof.
We only provide the proof for the first part, and the proof for the other parts is similar. For map (2), there is a fixed point with a multiplier +1, if(10)Φx,y,ϑ=x,yT,detJx,y−I2=0,where(11)I2=1001.
The center manifold atλ=r+d2 can be considered as(12)Mv=qv+m2v2+Ov3,M:ℝ⟶ℝ2,m2=m21,m22T.
The restriction of map (2) to the center manifold has the following form:(13)v⟶v+aLPv2+Ov3,where Jq=q,JTp=p,<p,q>=1, and(14)q=−d2d11,p=01.
By invariance property of the center manifold, we conclude that(15)aLP=12p,Bq,q=−hd1r+d2A.
Note thatE0 is always the fixed point and it will not be destroyed and aLP≠0. It is concluded that the fixed point E0 undergoes a transcritical bifurcation.Theorem 2.
(1)
There is a nondegenerate flip bifurcation ofE0 at λ=hr+hd2−2/h, when hd1≠2(2)
There is a nondegenerate flip bifurcation ofE0 at h=2/r+d2−λ, when hd1≠2(3)
There is a nondegenerate flip bifurcation ofE0 at d2=2+hλ−hr/h, when hd1≠2(4)
There is a nondegenerate flip bifurcation ofE0 at r=2+hλ−hd2/h, when hd1≠2Proof.
We only provide the proof for the first part, and the proof for the other parts is similar. A fixed point of model (2) with the multiplier −1 satisfies the following algebraic system:(16)Φx,y,ϑ=x,yT,detJx,y+I2=0.
The center manifold atλ=hr+hd2−2/h can be considered as(17)Mv=qv+m2v2+m3v3+Ov4,M:ℝ⟶ℝ2,mi=mi1,mi2T,i=2,3.
The restriction of map (2) to the center manifold has the form(18)ζ⟶−ζ+bPDζ3+Oζ4,where Jq=−q,JTp=−p,p,q=1, and(19)q=−hd2−2hd1−21,p=01.
As regards the invariance property of the center manifold, the critical normal form coefficients of the flip bifurcation can be presented as follows:(20)bPD=16p,Cq,q,q+3Bq,In−J−1Bq,q=d12hr+hd2−2A2hd1−2hd1−2hr+hd2−2−hd2−d1.
The flip bifurcation is nondegenerate providedbPD≠0. If bPD>0, the flip bifurcation is supercritical and a stable double-period cycle bifurcates from E0. If bPD<0, the flip bifurcation is subcritical and the double-period cycle is unstable.


### 3.2. Codimension One Bifurcation atE1

Theorem 3.
(1)
There is a nondegenerate flip bifurcation ofE1 at λ=2hd1−hr+d22+hd1−2hd2−4−G1/2hhd2−2, where G1=hd1−26h2d22+4hd2hr−2−22+hr2+hd1hr−2+hd22(2)
There is a nondegenerate flip bifurcation ofE1 at h=λλ+d1−r−d2−G2/λ−r−d2d2λ+d1−r−d2+rd1, where G2=λλd1+λ−4d2λ−r−d22+2d1λ−r−d2λ−2r−2d2(3)
There is a nondegenerate flip bifurcation ofE1 at r=2λ+hd1λ+2hd2d2−λ−d1−G3/2hd1−d2, where G3=λhd1−28d2−8d1−2λ+λhd1(4)
There is a nondegenerate flip bifurcation ofE1 at d1=hλ−r−d2hd2r+d2−λ+2λ−4λ/hr+d2λ−r−d2−2λhProof.
We only provide the proof for the first part, and the proof for the other parts is similar.
It is clear that at the pointE1, map (2) has a fixed point with a multiplier −1, or on the other hand, this point convinces the following equation:(21)Φx,y,ϑ−x,yT=0,detJx,y+I2=0.
The corresponding characteristic equation atE1 is given by(22)ρ2−trJρ+detJ=0,where(23)trJ=2+hr+d2−λ−d1,detJ=1+hr+d2−λ−d1+h2λ−r−d2d2λ+d1−r−d2+rd1λ.
The next proof is similar to the proof of Theorem2.
To avoid the complexity of the computation, we set(24)A=8,d1=0.08,d2=0.1,r=0.1,h=3.76.
Then, we have(25)q=−0.95951,p=01.
By using the invariance property of the center manifold, we can obtain the critical normal form coefficientbPD=0.0061. Since bPD>0, the bifurcation is supercritical and a stable double-period cycle bifurcates from E1.Theorem 4.
(1)
There is a nondegenerate Neimark–Sacker bifurcation ofE1 at h=λλ+d1−r−d2/λ−r−d2d2λ+d1−r−d2+rd1.(2)
There is a nondegenerate Neimark–Sacker bifurcation ofE1 at λ=d1−r+d2hd1−2hd2+1+H2/2hd2−1, where H2=r+d22+2d1r+d2hd2−hr−1+d12hd2+hr−12.(3)
There is a nondegenerate Neimark–Sacker bifurcation ofE1 at r=2hd2d2−d1−λ+hλd1+λ+H3/2hd1−d2, where H3=λλ+2hd12d2−λ+hd12−4+hλ.(4)
There is a nondegenerate Neimark–Sacker bifurcation ofE1 at d1=λ−r−d2λ−hd2λ−r−d2/hr+d2λ−r−d2−λ.Proof.
Only the first part is proved, and the proof of the other part is similar.
Map (2) has a fixed point with a pair of complex conjugate multipliers on the unit circle if(26)Φx,y,ϑ=x,yT,detJx,y=1.
Let us consider(27)Mv,v¯=∑1≤k+l1k+l!mklβvkv¯l,v∈ℂ,mkl∈ℂ,as the center manifold at h=λλ+d1−r−d2/λ−r−d2d2λ+d1−r−d2+rd1. Map (2) can be restricted to the center manifold which at the critical value is as follows:(28)ζ⟶eiθ0ζ+dNSζζ2+Oζ4,ζ∈ℂ,where dNS is a complex number. To avoid the complexity of the computation, we let(29)A=6,d1=0.145,d2=0.42,r=0.12,λ=0.75.
Thus,E1 has a pair of a simple multiplier.(30)ρ1,2=−0.3516±0.9362i.
This applies to the nonresonance conditions. Note that(31)q=−0.4473+2.0907i1,p=0.4266−0.1758i−0.1768−0.9705i,have been used, where(32)Jq=eiθ0q,JTp=e−iθ0p,p,q=1.
And,dNS is a complex number given by(33)dNS=12p,Cq,q,q¯+2Bq,I2−J−1Bq,q¯+Bq¯,e2iθI2−J−1Bq,q.
And,cNS=ℜe−iθdNS.
After long but basic computations, we can achievedNS=−0.1315−0.2105i; therefore, cNS=−0.1508. Since cNS<0, the Neimark–Sacker bifurcation is supercritical and the closed curve is stable.


## 3.1. Codimension One Bifurcation atE0

Theorem 1.
(1)
There is a transcritical bifurcation ofE0 at λ=r+d2 when hd1≠2(2)
There is a transcritical bifurcation ofE0 at r=λ−d2 when hd1≠2(3)
There is a transcritical bifurcation ofE0 at d2=λ−r when hd1≠2Proof.
We only provide the proof for the first part, and the proof for the other parts is similar. For map (2), there is a fixed point with a multiplier +1, if(10)Φx,y,ϑ=x,yT,detJx,y−I2=0,where(11)I2=1001.
The center manifold atλ=r+d2 can be considered as(12)Mv=qv+m2v2+Ov3,M:ℝ⟶ℝ2,m2=m21,m22T.
The restriction of map (2) to the center manifold has the following form:(13)v⟶v+aLPv2+Ov3,where Jq=q,JTp=p,<p,q>=1, and(14)q=−d2d11,p=01.
By invariance property of the center manifold, we conclude that(15)aLP=12p,Bq,q=−hd1r+d2A.
Note thatE0 is always the fixed point and it will not be destroyed and aLP≠0. It is concluded that the fixed point E0 undergoes a transcritical bifurcation.Theorem 2.
(1)
There is a nondegenerate flip bifurcation ofE0 at λ=hr+hd2−2/h, when hd1≠2(2)
There is a nondegenerate flip bifurcation ofE0 at h=2/r+d2−λ, when hd1≠2(3)
There is a nondegenerate flip bifurcation ofE0 at d2=2+hλ−hr/h, when hd1≠2(4)
There is a nondegenerate flip bifurcation ofE0 at r=2+hλ−hd2/h, when hd1≠2Proof.
We only provide the proof for the first part, and the proof for the other parts is similar. A fixed point of model (2) with the multiplier −1 satisfies the following algebraic system:(16)Φx,y,ϑ=x,yT,detJx,y+I2=0.
The center manifold atλ=hr+hd2−2/h can be considered as(17)Mv=qv+m2v2+m3v3+Ov4,M:ℝ⟶ℝ2,mi=mi1,mi2T,i=2,3.
The restriction of map (2) to the center manifold has the form(18)ζ⟶−ζ+bPDζ3+Oζ4,where Jq=−q,JTp=−p,p,q=1, and(19)q=−hd2−2hd1−21,p=01.
As regards the invariance property of the center manifold, the critical normal form coefficients of the flip bifurcation can be presented as follows:(20)bPD=16p,Cq,q,q+3Bq,In−J−1Bq,q=d12hr+hd2−2A2hd1−2hd1−2hr+hd2−2−hd2−d1.
The flip bifurcation is nondegenerate providedbPD≠0. If bPD>0, the flip bifurcation is supercritical and a stable double-period cycle bifurcates from E0. If bPD<0, the flip bifurcation is subcritical and the double-period cycle is unstable.


## 3.2. Codimension One Bifurcation atE1

Theorem 3.
(1)
There is a nondegenerate flip bifurcation ofE1 at λ=2hd1−hr+d22+hd1−2hd2−4−G1/2hhd2−2, where G1=hd1−26h2d22+4hd2hr−2−22+hr2+hd1hr−2+hd22(2)
There is a nondegenerate flip bifurcation ofE1 at h=λλ+d1−r−d2−G2/λ−r−d2d2λ+d1−r−d2+rd1, where G2=λλd1+λ−4d2λ−r−d22+2d1λ−r−d2λ−2r−2d2(3)
There is a nondegenerate flip bifurcation ofE1 at r=2λ+hd1λ+2hd2d2−λ−d1−G3/2hd1−d2, where G3=λhd1−28d2−8d1−2λ+λhd1(4)
There is a nondegenerate flip bifurcation ofE1 at d1=hλ−r−d2hd2r+d2−λ+2λ−4λ/hr+d2λ−r−d2−2λhProof.
We only provide the proof for the first part, and the proof for the other parts is similar.
It is clear that at the pointE1, map (2) has a fixed point with a multiplier −1, or on the other hand, this point convinces the following equation:(21)Φx,y,ϑ−x,yT=0,detJx,y+I2=0.
The corresponding characteristic equation atE1 is given by(22)ρ2−trJρ+detJ=0,where(23)trJ=2+hr+d2−λ−d1,detJ=1+hr+d2−λ−d1+h2λ−r−d2d2λ+d1−r−d2+rd1λ.
The next proof is similar to the proof of Theorem2.
To avoid the complexity of the computation, we set(24)A=8,d1=0.08,d2=0.1,r=0.1,h=3.76.
Then, we have(25)q=−0.95951,p=01.
By using the invariance property of the center manifold, we can obtain the critical normal form coefficientbPD=0.0061. Since bPD>0, the bifurcation is supercritical and a stable double-period cycle bifurcates from E1.Theorem 4.
(1)
There is a nondegenerate Neimark–Sacker bifurcation ofE1 at h=λλ+d1−r−d2/λ−r−d2d2λ+d1−r−d2+rd1.(2)
There is a nondegenerate Neimark–Sacker bifurcation ofE1 at λ=d1−r+d2hd1−2hd2+1+H2/2hd2−1, where H2=r+d22+2d1r+d2hd2−hr−1+d12hd2+hr−12.(3)
There is a nondegenerate Neimark–Sacker bifurcation ofE1 at r=2hd2d2−d1−λ+hλd1+λ+H3/2hd1−d2, where H3=λλ+2hd12d2−λ+hd12−4+hλ.(4)
There is a nondegenerate Neimark–Sacker bifurcation ofE1 at d1=λ−r−d2λ−hd2λ−r−d2/hr+d2λ−r−d2−λ.Proof.
Only the first part is proved, and the proof of the other part is similar.
Map (2) has a fixed point with a pair of complex conjugate multipliers on the unit circle if(26)Φx,y,ϑ=x,yT,detJx,y=1.
Let us consider(27)Mv,v¯=∑1≤k+l1k+l!mklβvkv¯l,v∈ℂ,mkl∈ℂ,as the center manifold at h=λλ+d1−r−d2/λ−r−d2d2λ+d1−r−d2+rd1. Map (2) can be restricted to the center manifold which at the critical value is as follows:(28)ζ⟶eiθ0ζ+dNSζζ2+Oζ4,ζ∈ℂ,where dNS is a complex number. To avoid the complexity of the computation, we let(29)A=6,d1=0.145,d2=0.42,r=0.12,λ=0.75.
Thus,E1 has a pair of a simple multiplier.(30)ρ1,2=−0.3516±0.9362i.
This applies to the nonresonance conditions. Note that(31)q=−0.4473+2.0907i1,p=0.4266−0.1758i−0.1768−0.9705i,have been used, where(32)Jq=eiθ0q,JTp=e−iθ0p,p,q=1.
And,dNS is a complex number given by(33)dNS=12p,Cq,q,q¯+2Bq,I2−J−1Bq,q¯+Bq¯,e2iθI2−J−1Bq,q.
And,cNS=ℜe−iθdNS.
After long but basic computations, we can achievedNS=−0.1315−0.2105i; therefore, cNS=−0.1508. Since cNS<0, the Neimark–Sacker bifurcation is supercritical and the closed curve is stable.


## 4. Codimension Two Bifurcation Analysis

In this section, we study the codimension two bifurcations of model (2). By using the Jacobian matrix and some strict mathematical calculations, we find that E0 does not come into being any codimension two bifurcations.ForE1, the following theorems give the codimension two bifurcation analysis.Theorem 5.
(1)
Ifd1=λ−r−d2r−d2+r+d22−rλ/r+d2−r+d22−rλ and h=2r+2d2−2r+d22−rλ/rλ−r−d2, there is a 1 : 2 strong resonance(2)
Ifd1=2hr−hd2+1+2hr+d2hd2−2−1/rh2 and λ=hr+d22+hr−2−2hr+d2hd2−2−1/rh2, there is a 1 : 2 strong resonance(3)
Ifd1=22hr−hλ+8+2hλhλ−hr−4/h4+hr and d2=2hλ−2r+rh2λ−r+2hλhλ−hr−4/h4+hr, there is a 1 : 2 strong resonanceProof.
Only the first part is proved, and the proof of the other parts is similar.
Map (2) experiences a 1 : 2 strong resonance if it has a fixed point with two multipliers −1. Then, we have(34)Φx,y,ϑ=x,yT,detJx,y=1,trJx,y=−2.
The exact solution is as follows:(35)x=Ar+d2r+d2−r+d22−rλrλ−r−d2r+d2+r+d22−rλ,y=Ar+d2−r+d22−rλrr+d2+r+d22−rλ,h=2r+2d2−2r+d22−rλrλ−r−d2,d1=λ−r−d2r−d2+r+d22−rλr+d2−r+d22−rλ.
Let us consider(36)Mv1,v2=v1q0+v2q1+∑2≤j+k<31j!k!mjkv1jv2k,as the center manifold at E1. Map (2) can be restricted to the center manifold at the critical value parameters of handd1 as follows:(37)ξ1ξ2⟶−ξ1+ξ2−ξ2+CR2ξ13+DR2ξ12ξ2,in which(38)Jq0=−q0,Jq1=−q1+q0,JTp0=−p0,JTp1=−p1+p0,p0,q1=p1,q0=1,p0,q0=p1,q1=0.
We have(39)q0=−rλ+r+d2r+d2−r+d22−rλλ−r−d2r+d2−r+d22−rλ1,q1=rλλ−r−d2r+d2−r+d22−rλ0,p0=λ−r−d2r+d2−r+d22−rλrλrλ−r+d2r+d2−r+d22−rrλ,p1=01,CR2=16p0,Cq0,q0,q0+3Bq0,h20,DR2=12p0,2Bq0,h11+Bq1,h20+Cq0,q0,q1+12p1,3Bq0,h20+Cq0,q0,q0,where h20,h11,andh02 satisfy(40)J−I2h20=−Bq0,q0,J−I2h11=−Bq0,q1−h20.
To avoid the complexity of the computation, we letA=6,r=0.12,d2=0.42,andλ=0.75. After long but basic computations, we obtain the explicit solutions C1R2=4CR2=0.3012≠0 and D1R2=−2DR2−6CR2=−2.0864≠0, which are the nondegenerate conditions of 1 : 2 resonance. The sign of C1R2 determines the criticality of this bifurcation. The bifurcation scenario is indicated by the coefficient DR2.Theorem 6.
(1)
Ifd1=λ−r−d23r−3d2+9r+d22−12rλ/3r+3d2−9r+d22−12rλ and h=3r+3d2−9r+d22−12rλ/2rλ−r−d2, there is a 1 : 3 strong resonance(2)
Ifd1=3hd2+hr+1−3hr+d23hd2−hr−6+9/2rh2 and λ=hr+d23+2rh−3+3hr+d23hd2−hr−6+9/2rh2, there is a 1 : 3 strong resonance(3)
Ifd1=32hr−hλ+6+3hλ3hλ−4hr−12/2h3+hr and d2=3hλ−2r+2rh2λ−r+3hλ3hλ−4hr−12/2h3+hr, there is a 1 : 3 strong resonanceProof.
Only the first part is proved, and the proof of the other parts is similar.
Map (2) undergoes a 1 : 3 strong resonance if it has a fixed point with a pair of conjugate multipliers e±i2π/3. Then,(41)Φx,y,ϑ=x,yT,detJx,y=1,trJx,y=−1.
The exact solution is as follows:(42)x=Ar+d23r+3d2−9r+d22−12rλrλ−r−d23r+3d2+9r+d22−12rλ,y=A3r+3d2−9r+d22−12rλr3r+3d2+9r+d22−12rλ,h=3r+3d2−9r+d22−12rλ2rλ−r−d2,d1=λ−r−d23r−3d2+9r+d22−12rλ3r+3d2−9r+d22−12rλ.
The restriction of model (2) to two-dimensional center manifolds is expressed as follows:(43)Mv,v¯=vq+v¯q¯+∑2≤j+k<31j!k!mjkvjv¯k.
Let(44)q=−3−3irλ−r+d23r+3d2−9r+d22−12rλλ−r−d23r+3d2−9r+d22−12rλ1,p=3i−12λ−r−d23r+3d2−9r+d22−12rλ3−3irλ−r+d23r+3d2−9r+d22−12rλ1+3i2.
It is easy to check that(45)Jq=e2π/3iq,JTp=e−2π/3ip,p,q=1.
Model (2) can be transformed into the following normal form:(46)ξ⟶e2π/3iξ+BR3ξ¯2+CR3ξξ2+Oξ4,where BR3andCR3 can be determined by(47)BR3=12p,Bq¯,q¯,CR3=12p,Cq,q,q¯+2Bq,h11+Bq¯,h20.
We have(48)J−e4π/3iI2h20=2B¯R3q¯−Bq,q,J−I2h11=−Bq,q¯.
To avoid the complexity of the computation, we letA=6,r=0.12,d2=0.42,andλ=0.75. Then, we have BR3=0.0918−0.3755i,andCR3=0.4313+0.7614i, and the stability of the bifurcating closed invariant curve is determined by ℜ3e4π/3iCR3−BR32.Theorem 7.
(1)
Ifd1=λ−r−d2r−d2+r+d22−2rλ/r+d2−r+d22−2rλ and h=r+d2−r+d22−2rλ/rλ−r−d2, there is a 1 : 4 strong resonance(2)
Ifd1=1−2hd2+2hr+d22+hrhr+hd2−1+1/rh2 and λ=hr+d22+rh−1−2hr+d22+hrhr+hd2−1+1/rh2, there is a 1 : 4 strong resonance(3)
Ifd1=2hr−hλ+4+hλhλ−2hr−4/h2+hr and d2=hλ−2r+rh2λ−r+hλhλ−2hr−4/h2+hr, there is a 1 : 4 strong resonanceProof.
Only the first part is proved, and the proof of the other parts is similar.
Map (2) experiences a 1 : 4 strong resonance if it has a fixed point with a pair of conjugate multipliers ±i. Then,(49)Φx,y,ϑ=x,yT,detJx,y=1,trJx,y=0.
The exact solution is as follows:(50)x=Ar+d2r+d2−r+d22−2rλrλ−r−d2r+d2+r+d22−2rλ,y=Ar+d2−r+d22−2rλrr+d2+r+d22−2rλ,h=r+d2−r+d22−2rλrλ−r−d2,d1=λ−r−d2r−d2+r+d22−2rλr+d2−r+d22−2rλ.
The center manifold atE1 can be considered as follows:(51)Mv,v¯=vq+v¯q¯+∑2≤j+k<31j!k!mjkvjv¯k.
Let(52)q=−rλ1−i−r+d2r+d2−r+d22−2rλλ−r−d2r+d2−r+d22−2rλ1,p=i−12λ−r−d2r+d2−r+d22−2rλrλ1−i−r+d2r+d2−r+d22−2rλ1+i2.
It is easy to check that(53)Jq=iq,JTp=−ip,p,q=1.
Model (2) can be transformed into the following normal form:(54)ξ⟶iξ+CR4ξξ2+DR4ξ¯3+Oξ4,where CR4andDR4 can be achieved by the following equation:(55)CR4=12p,Cq,q,q¯+2Bq,I2−J−1Bq,q¯−Bq¯,I2+J−1Bq,q,DR4=16p,Cq¯,q¯,q¯−3Bq¯,I2+J−1Bq¯,q¯.
The bifurcation scenario near the 1 : 4 resonance point is determined by(56)C0=−iCR4DR4.
To avoid the complexity of the computation, we letA=6,r=0.12,d2=0.42,andλ=0.75. Then, we have C0=0.0938−0.194i.


## 5. Numerical Simulations

In this section, we will give bifurcation diagrams, phase portraits, and Lyapunov exponents to illustrate the above theoretical analyses and find new interesting complex dynamical behaviors by using numerical simulations. The bifurcation parameters are considered in the following two cases.Case 1.
Varyh in range h∈3.6,5.2, and fix A=8,d1=0.08,d2=0.1,r=0.1,andλ=0.75. By calculation, we find that at a unique endemic equilibrium E122.5352,61.9728, a flip bifurcation occurs at h=3.7583.
The flip bifurcation mentioned in Theorem3 is depicted in Figure 1(a), from which we see that stability of fixed point E1 happens for h<3.7583 and loses its stability at h=3.7583 and period-doubling phenomena lead to chaos for h>3.7583. The maximum Lyapunov exponent related to Figure 1(a) is computed and shown in Figure 1(b). We observe that the period 2,4,8,16 orbits occur for h∈3.7583,4.792 and chaotic sets for h∈4.792,5.2. The period 6 orbit occurs at h=4.932 and period 5 orbit occurs at h=5.104 within the window of the chaotic region. Here, the powerful Lyapunov exponent test is shown in Figure 1(b) to confirm the existence of periodic orbits and chaotic motions.Figure 1
(a) Flip bifurcation diagram of map (2). (b) The maximum Lyapunov exponent.
(a)(b)Case 2.
FixA=6,d1=0.145,d2=0.42,r=0.12,andλ=0.75 and choose h∈7,8.6. By calculation, we observe that a Neimark–Sacker (NS) bifurcation appears at a unique epidemic equilibrium E119.459459,7.567568 for h=7.6148. Figure 2 shows the correctness of Theorem 4.
The bifurcation diagrams shown in Figures2(a) and 2(b) demonstrate that stability of E1 happens for h<7.6148 and there is a loss in its stability at h=7.6148 and an attracting invariant curve appears if h>7.6148. We dispose the maximum Lyapunov exponent in Figure 2(c) relating bifurcation diagrams, which confirms the existences of chaos and period windows as parameter h varies. When h∈8.202,8.6, the sign of the maximum Lyapunov exponent is greater than zero, confirming the presence of chaos. The status of stable, periodic, or chaotic dynamics is compatible with the sign in Figure 2(a). As we show, all our analytical predications have excellent arguments with numerical results.
The dynamical complexity of model (2) can be observed when more parameters vary in Figure 3. The bifurcation diagrams of system (2) for controlling parameters h∈7,8.3 and d1∈0.13,0.18 and fixing remaining parameters as in Case 2 are shown in Figures 3(a) and 3(b). It is easy to find values of control parameters for which the dynamics of system (2) are in the status of quasiperiodic, periodic, or chaotic. This shows that the parameter d1 may play an important role in chaotic dynamics in this system.
By the mathematical results in Theorem6, we know that there exists 1 : 3 resonance when h=7.4861andd1=0.1907 with x,y=16.945541,6.589933, whose parameters are kept fixed at A=6,d2=0.42,r=0.12,andλ=0.75. To prove it, we present a global view of the d1×h parameter space in Figure 4(a), for 0<d1<0.28 and 7<h<8.6. Figure 4(a) is a diagram obtained by discretizing the parameter interval in a grid of 1000×1000 points equally spaced [23, 24]. For each point h,d1, an orbit of initial condition x0,y0=16.95,6.59 converges to a chaotic, or a quasiperiodic, or a periodic attractor, or an attractor at infinity (diverges). In Figure 4(a), the period −1 solution is shown by orange, the period −29 is plotted in light gray, periods 2−28 are codified in other different colors, and period motions are indicated by numbers. Moreover, quasiperiodic regions are indicated by the letter Q. Divergence regions are indicated by the letter D. Chaotic regions are indicated by the letter C. From Figure 4(a), we can also find 1 : 3 resonance when choosing a dark cyan region near the neighbourhood of h=7.4861andd1=0.1907. This is entirely consistent with the results we obtain analytically.
Figure4(b) shows the magnification of a region with quasiperiodic behaviors in Figure 4(a) with 7.5<h<8.6 and 0.09<d1<0.16, which are periodic regions similar to Arnold tongues [25, 26]. For points on the white-orange borderline in Figure 4(b), on which Neimark–Sacker bifurcation (NSb) occurs, period −1 orbits disappear, to give rise to invariant closed curves, which are the characteristics of a quasiperiodic motion. Therefore, quasiperiodicity is related to parameters near the white-orange border, on the white side, while chaotic solutions are corresponding to parameters far away from this same border, also on the white side. As is known [25, 26], the periodic structures embedded in the quasiperiodic regions possess self-similarity fractal features. In addition, these periodic structures are numbered as 7, 10, 13, 16, 19, 22, and 25. As there exists 1 : 3 resonance in this region, the difference between the adjacent period is 3. Due to the period-doubling bifurcation, we can see Arnold tongue sequence 22, 14, 20, 26. This is the main message of this paper: period-adding cascades exist in profusion in the SIS model.
To our surprise, by checking the rotation number of these Arnold tongues in Figure4(b), it is observed that the mode-locking order is organized according to the broken Farey tree, which is organized according to a well-defined rule in [26]. Between each generic tongue of rotation numbers p/q and r/s, there is a tongue with a period equal to q+s and rotation number p+r/q+s. For instance, the tongue with the rotation number 5/18 is located between the tongues with the rotation numbers 3/11 and 2/7. Also are shown in Figure 4(b) the tongues 25 between 7 and 18, 24 between 7 and 17, 27 between 17 and 10, and 23 between 10 and 13. The sequence formed by the corresponding rotation numbers, namely, 5/18, 7/25, 2/7, 7/24, 5/17, 8/27, 3/10, 7/23, and 4/13, is only a small and incomplete part of a Farey tree, the so-called broken Farey tree.
In the sequel, we will concentrate on Figure4(b), to investigate details on the dynamics of model (2). Figure 2 shows the bifurcation diagram and the maximum Lyapunov exponent diagram when d1=0.145. By increasing h from 7 in Figure 2(a), the system stays in a period −1 state until h=7.6148, being characterized by the largest Lyapunov exponent less than zero. At h=7.6148, the largest Lyapunov exponent reaches the zero value and a quasiperiodic region is attained, remaining until h=7.977. The quasiperiodic motion is, therefore, originated from a period −1 orbit and is a consequence of an NSb that occurs at h=7.6148. As a result of this NSb, one closed invariant is created, whose typical example can be seen in Figure 5(a), which shows a phase portrait for h=7.78andd1=0.145. There are mode-locking states (Arnold tongues) embedded in the quasiperiodic regions, for which the largest Lyapunov exponent is less than zero. By increasing h from 7.6148, we cross the quasiperiodic region and penetrate a period −13 tongue at h=7.977. Period-doubling bifurcation happens at h=8.148, resulting in period −26 orbits. Examples of these periodic orbits can be seen in Figures 5(b) and 5(c). By increasing h more, finally, a chaotic state for which the largest Lyapunov exponent is greater than zero is created at h=8.173. A typical chaotic attractor appears in Figure 5(d), for h=8.205andd1=0.145.
From Theorem6, we know that there also exists 1 : 3 resonance when λ=0.745andd1=0.189 and d2=0.426andd1=0.191, whose parameters are fixed as A=6,d2=0.42,r=0.12,andh=7.6, and A=6,λ=0.75,r=0.12,andh=7.6, respectively. To prove it, we present the global view of the d1×λ and d1×d2 parameter planes in Figure 6, where the color representations are as the same as in Figure 4. Similar to Figure 4, there also exist periodic structures which are numbered as 7, 10, 13, 16, 19, and 22. As there exists 1 : 3 resonance in the dark cyan region, the difference between the adjacent period is 3. In addition, we can also find periodic motion, quasiperiodic motion, and chaotic motion appear alternately. And, the complex mode-locking structures can also be seen in Figure 6.Figure 2
(a) Neimark–Sacker bifurcation inh,x space. (b) Neimark–Sacker bifurcation in h,y space. (c) Maximum Lyapunov exponent related to a. (d) Magnification of (a) for 7.965<h<8.185.
(a)(b)(c)(d)Figure 3
Model (2) by varying the dynamical parameters handd1 and setting the parameter values A=6,r=0.12,d2=0.42,andλ=0.75 shows (a) the bifurcation diagram in h−d1−x space and (b) the bifurcation diagram in h−d1−y space.
(a)(b)Figure 4
(a) Global view of thed1,h parameter space, showing different stability regions. (b) Magnification of the rectangle in (a) is located in the range of 7.5<h<8.6 and 0.09<d1<0.16. In each diagram, the number is related to the period.
(a)(b)Figure 5
Phase portraits for various values ofh. (a) h=7.78, (b) h=8.055, (c) h=8.16, and (d) h=8.205.
(a)(b)(c)(d)Figure 6
(a) Thed1×λ parameter-space diagram for A=6,d2=0.42,r=0.12,andh=7.6 with the initial condition x0,y0=17.22,6.54. (b) The d1×d2 parameter-space diagram for A=6,λ=0.75,r=0.12,andh=7.6 with initial condition x0,y0=17.135,6.402.
(a)(b)


## 6. Conclusions

In summary, we give a detailed analysis for the codimension one bifurcation and codimension two bifurcation of a discrete-time SIS epidemic model. Although it is a simple discrete epidemic model, we find that there are many complex dynamical behaviors. The theoretical analysis demonstrates that model (2) undergoes codimension one bifurcation (transcritical, filp, Neimark–Sacker bifurcation) and codimension two bifurcation (1 : 2, 1 : 3, 1 : 4 strong resonances). For each bifurcation, the normal forms for maps and their approximations to the flows of the corresponding differential equations are obtained to determine the bifurcation diagrams. A good correlation between the experimental and theoretical results is found. Furthermore, we numerically report a visual view of the mode-locking structures to the discrete model in two-dimensional parameter-space diagrams. And, we also show Arnold tongue structures embedded in quasiperiodic regions regardless of any combination of parameters, respectively, organized themselves in the period-adding bifurcation cascades. The results obtained in this paper are interesting and useful and may be applied to a variety of epidemic models.



---

*Source: 2233452-2022-04-23.xml*