# Bifurcation and Chaotic Behavior of a Discrete-Time SIS Model

Authors: Junhong Li; Ning Cui
Journal: Discrete Dynamics in Nature and Society
 (2013)

---

## Abstract

The discrete-time epidemic model is investigated, which is obtained using the Euler method. It is verified that there exist some dynamical behaviors in this model, such as transcritical bifurcation, flip bifurcation, Hopf bifurcation, and chaos. The numerical simulations, including bifurcation diagrams and computation of Lyapunov exponents, not only show the consistence with the theoretical analysis but also exhibit the rich and complex dynamical behaviors.


---

## Body

## 1. Introduction

Epidemic models have been widely used in different forms for studying epidemiological processes such as the spread of HIV [1], SARS [2], and influenza [3]. It is well known that dynamical systems with simple dynamical behavior in the constant parameter case display very complex behaviors including chaos when they are periodically perturbed [4, 5]. The continuous-time epidemic models have been widely investigated in many articles (e.g., [6–10] and the references cited therein). In recent years, we have found that more attention is paid to the discrete-time epidemic models (see [11–15] and the references cited therein). The reasons are as follows: first, difference models are more realistic than continuous differential ones because the epidemic statistics are compiled from given time intervals and are discontinuous. Second, the discrete-time models can provide natural simulators for the continuous cases. One can thus not only study the behaviors of the continuous-time model with good accuracy, but also assess the effect of lager time steps. At last, the use of discrete -time models makes it possible to use the entire arsenal of methods recently developed for the study of mappings and lattice equations, either from the integrability and/or chaos points of view.On the other hand, daily treatments are frequently done for some infections, such as the group of those being responsible for the common cold, which do not confer any long lasting immunity. Such infections do not have a recovered state and individuals become susceptible again after infection. For such reasons, according to [16], we firstly consider the SIS epidemic model with nonlinear incidence rate:
(1)dSdt=bS(1-SK)-βS2I+γI,dIdt=βS2I-dI-γI,
where S denotes the susceptible population, I is the infected population, and b is the intrinsic birth rate constant. K, β are the carrying capacity and the infection rate, respectively. γ≥0 is the recovery rate constant (1/γ is the average infective time).Let(2)S=d+γβx,I=d+γβy,dτ=(d+γ)dt.
We obtain the following system analogous to (1):
(3)dxdτ=a0x-a1x2-x2y+a2y,dydτ=x2y-y,
where
(4)a0=bd+γ,a1=bKβ(d+γ),a2=γd+γ.
Applying Euler scheme to (3), we obtain the following equation:
(5)xn+1=(a0+1)xn-a1xn2-xn2yn+a2yn,yn+1=xn2yn.This paper is organized as follows. In Section2, we give sufficient conditions of existence for transcritical bifurcation, flip bifurcation, and Hopf bifurcation. In Section 3, a series of numerical simulations show that there are bifurcation and chaos in the discrete-time epidemic model. Finally, we give remarks to conclude this paper in Section 4.


## 2. Bifurcations

It is easy to visualize that system (5) has three fixed points P0(0,0), P1(a0/a1,0), and P2(1,(a0-a1)/(1-a2)) when a1, a0 are fixed. We can see that the fixed point P0(0,0) is a saddle. In the following paper, we focus on investigating the bifurcations of P1, P2.Theorem 1.
Ifa0=a1 and a0≠2, (5) undergo a transcritical bifurcation at P1. Furthermore, the system has three fixed points when a0>a1 and has two fixed points when a0≤a1.Proof.
The Jacobian matrix of (5) at P1 takes the form
(6)J(P1)=(1-a0a2-a02a120a02a12).J(P1) has eigenvalues λ1=1-a0, λ2=a02/a12. And a0≠2 implies that |λ1|≠1.
Let(7)ξn=xn-a0a1,ηn=yn,cn=a0-a1,
equation (5) becomes
(8)ξn+1=(1-a1)ξn+(a2-1)ηn+F1(ξn,ηn,cn),ηn+1=ηn+F2(ξn,ηn,cn),cn+1=cn,
where
(9)F1(ξn,ηn,cn)=-(a1+ηn)ξn2-2ξnηncna1-(2ηn+cn)ξn-ηncn2a12-2cnηna1,F2(ξn,ηn,cn)=2ξnηn+2cnηna1+ξn2ηn+ηncn2a12+2ξnηncna1.
By the following transformation:
(10)(ξnηncn)=(1101a1a2-10001)(pnqnσn),
equation (8) becomes
(11)(pn+1qn+1σn+1)=(1-a100010001)(pnqnσn)+(f1(pn,qn,σn)f2(pn,qn,σn)0),
where
(12)f1(pn,qn,σn)=(a1+a2-11-a2qn-a1)pn2+[2(a1+a2-1)(σn+1)(qn+qn2)a1(1-a2)-σn-2a1qn]pn+2(a1+a2-1)qn3a1(1-a2)+[2(a1+a2-1)(σn+a1)a1(1-a2)-a1]qn2+[(a1a2+2a2-2)σna1(1-a2)+2(a1+a2-1)σn2a13(1-a2)]qn,f2(pn,qn,σn)=qn3+2(1+pn+σna1)qn2f2(pn,qn,σn)=+2pnqn+2qnσna1+(a1pn+σn)2qna12.
Then, we can consider
(13)pn=g(qn,σn)=ε1qn2+ε2qnσn+ε3σn2+o((|pn|+|qn|)3),
which must satisfy
(14)g(qn+1,σn+1)=g(qn+f2(g(qn,σn),qn,σn),σn+1)=(1-a1)g(qn,σn)+f1(g(qn,σn),qn,σn).
Thus, we have
(15)ε1=a1a2+2a2-2a1(1-a2),ε2=a1a2+2a2-2a12(1-a2),ε3=0.
And (8) is restricted to the center manifold, which is given by
(16)f:qn+1=2(a1a2+2a2-2)a1(1-a2)(qn3+qn2σna1)+qn3+2qn2+qn+2qn2σna1+2qnσna1+o((|qn|+|σn|)4).
Since
(17)f(0,σn)=0,∂f∂q|(0,0)=1,∂2f∂q2|(0,0)=4,∂2f∂q∂σ|(0,0)=-2a1≠0,
system (5) undergoes a transcritical bifurcation at P1. This proves the theorem.Theorem 2.
Equation (5) undergoes a flip bifurcation at P1 when a0=2, a1≠a0. Furthermore, the stable periodic-2 point bifurcates from this fixed point.Proof.
Let(18)ξn=xn-a0a1,ρn=yn,cn=a0-2,
system (5) becomes
(19)ξn+1=-ξn+(a2-4a12)ρn+F1(ξn,cn,ρn),cn+1=-cn,ρn+1=4a12ρn+F2(ξn,cn,ρn),
where
(20)F1(ξn,cn,ρn)=-(a1+ρn)ξn2-(cn+2(cn+2)ρna1)ξn-4cn+cn2a12ρn,F2(ξn,cn,ρn)=ρnξn2+4cn+cn2a12ρn+2(cn+2)ρna1ξn.
By the following transformation:
(21)(ξncnρn)=(10101000(a12+4)(a2a12-4))(ϕnunvn),
equation (19) becomes
(22)(ϕn+1un+1vn+1)=(-1000-10004a12)(ϕnunvn)+(f1(ϕn,un,vn)0f2(ϕn,un,vn)),
where
(23)f1(ϕn,un,vn)=(a12+a2a124-a2a12vn-a1)ϕn2+[2(1+a2)a12vn24-a2a12+2a1(un-2+a2a12+a2un+2a2)vn4-a2a12-un]ϕn+a12+a2a124-a2a12vn3+a1(a2a12+2a2un+4a2+2un)4-a2a12vn2+(1+a2)un2vn4-a2a12+a2(4+a12)unvn4-a2a12,f2(ϕn,un,vn)=vn3+2(ϕn+un+4a1)vn2f2(ϕn,un,vn)=+(a1ϕn+un)(a1ϕn+un+4)vna12.
Then, we can consider vn=g(ϕn,un)=ε1ϕn2+ε2ϕnun+ε3un2+o((|ϕn|+|un|)3), which must satisfy
(24)g(-ϕn+f1(ϕn,g(ϕn,un),vn),un+1)=4a12g(ϕn,un)+g3(ϕn,un)+2(ϕn+un+4a1)g2(ϕn,un)+(a1ϕn+un)(a1ϕn+un+4)g(ϕn,un)a12.
Thus, we have εk=0, k=1,2,3. We obtain the center manifold as follows:
(25)vn=g(ϕn,un)=ε4ϕn2un+ε5ϕnun2+ε6un3+o((|ϕn|+|un|)4).
And (19) is restricted to the center manifold, which is given by
(26)f:ϕn+1=-ϕn+f1(ϕn,g(ϕn,un),vn).
Direct calculations show that
(27)(∂f∂σ∂2f∂p2+2∂2f∂p∂σ)|(0,0)=-2,(12(∂2f∂p2)2+13(∂3f∂p3))|(0,0)=a122.
Hence, system (5) undergoes a flip bifurcation at P1. This completes the proof.The positive fixed point is so important to the biological system that people usually are very interested in it. We will next pay attention to the only positive fixed pointP2 of (5).Theorem 3.
Equation (5) undergoes a transcritical bifurcation at P2 when a1=a0, a1≠2, and Δ>0, where
(28)Δ=(1+a2)2a02(1-a2)2+8(a1-a0)+4a12a22-4a2(a2+1)(1-a2)2.Theorem 4.
IfΔ>0, a1≠a0, and a0=(2-2a2-a1+3a1a2)/2a2, (5) undergoes a flip bifurcation at P2.Since the analysis is similar to the case atP1, the previously mentioned proofs are omitted.We next give the condition of existence of Hopf bifurcation by using the hopf bifurcation theorem in [17].The characteristic equation of the Jacobian matrixJ(P2) can be written as λ2-t2λ+d2=0, where
(29)t2=2+a0-2a1-2(a0-a1)1-a2,d2=t2-1+2a0-2a1.The eigenvalues of the Jacobian matrix of (5) at P2 are λ1,2=(t2±t22-4d2)/2. The eigenvalues λ1,2 are complex conjugates for Δ<0. We translate the fixed point P2(1,(a0-a1)/(1-a2)) to the origin by xn=φn-1, yn=δn-((a0-a1)/(1-a2)), and the system (5) becomes
(30)φn+1=(t2-1)φn+(a2-1)δn-(δn+a0-a1a21-a2)φn2-2φnδn,δn+1=2(a0-a1)1-a2φn+δn+2φnδn+(δn+a0-a11-a2)φn2.
The eigenvalues of the matrix associated with the linearized map (30) at fixed point (0,0) are complex conjugates which are written as λ,λ-=(t2±i4d2-t22)/2, where
(31)|λ|=1+3a0-4a1-2(a0-a1)1-a2.
Now assume that 3a2<1, and let a00=(2a1(1-2a2))/(1-3a2). Then we have
(32)d(|λ|)da0|a0=a00=1-3a22(1-a2)>0,|λ(a10)|=1,λ(a00)=1+a1a2-a11-3a2±ia1(1-a2)(2-a1+a1a2-6a2)1-3a2,λj≠1,j=1,2,3,4.
Let
(33)T=(a2-102-t22-4d2-t222),(φnδn)=T(ϕnσn).
The system (30) becomes
(34)(ϕn+1σn+1)=(t22-4d2-t2224d2-t222t22)(ϕnσn)+(f1(ϕn,σn)f2(ϕn,σn)),
where
(35)f1(ϕn,σn)=(a0+a0a22-a1a2)ϕn3+[(a2-1)4d2-t222σn+a2(2a0-a1-a1a2)a2-1]ϕn2+4d2-t22σnϕn,f2(ϕn,σn)=(2-t2)(2a2-t2)(a2-1)24d2-t22ϕn3+(a2-1)(a2-t22)ϕn2σn-(2-t2)(2-t2+a1a2-a0)4d2-t22ϕn2-2(a2-1)(2-t2-a0+a1)4d2-t22ϕn2+(2a2-t2)σnϕn.
Notice that (34) is exactly on the center manifold in the form, in which the coefficient l [18] is given by l=-Re[((1-2λ)λ-2/(1-λ))l11l20]-(1/2)|l11|2-|l02|2+Re(λ-l21), where
(36)l11=14[(f1ϕϕ+f1σσ)+(f2ϕϕ+f2σσ)i]=3a0a2-2a1a22-a04a2-4-144d2-t22×[(4-2t2)(2-t2+a1a2-a0)+(a2-1)(4-2t2-2a0+2a1)]i,l20=18[(f1ϕϕ-f1σσ+2f2ϕσ)+(f2ϕϕ-f2σσ-2f1ϕσ)i]=14(a2-a1a2-2+2a0)-184d2-t22×[{t22}(4-2t2)(2-t2+a1a2-a0)+(a2-1)×(4-2t2-2a0+2a1)-8d2+2t22]i,l02=18[(f1ϕϕ-f1σσ-2f2ϕσ)+(f2ϕϕ-f2σσ+2f1ϕσ)i]=18(2a0-2a0a2-2-4a2+3t2)+18(2a0-a1a2-2+t2+24d2-t22)i,l21=116[(f1ϕϕϕ+f1ϕσσ+f2ϕϕσ+f2σσσ)+(f2ϕϕϕ+f2ϕσσ-f1ϕϕσ-f1σσσ)i]=18(1-a2)(3+a2-2t2)+184d2-t22(6a2-3t2+3a2t2+8d2+t22)i.
From the previous analysis, we have Theorem 5.Theorem 5.
System (5) undergoes a Hopf bifurcation at fixed point P2, if l≠0, Δ<0 and
(37)a1≠2(1-a2)+2a2a11+a3.


## 3. Numerical Simulations

With development of scientific computation, computer becomes a powerful tool to study nonlinear systems, especially for a system without explicit solution. It not only is able to explore new complex dynamical behaviors, for example, periodic orbits and chaos in different regions [16, 19], but also can do rigorous analysis by combining modern dynamical systems theory and reliable computation, for example, topological horseshoes [20–23]. In this section, we use the bifurcation diagrams, Lyapunov exponents, and phase portraits to illustrate the previous analytic results and find new dynamics of the model (5) as the parameters vary. The attractors of (5) are also given by using the method in [20]. The bifurcation parameters are considered in the following three cases:(I)
varyinga0 in range 1.6≤a0≤2.2 and fixing a1=0.9, a3=0.01;(II)
varyinga1 in range 0.8≤a1≤1.2 and fixing a0=2, a3=0.01;(III)
varyinga2 in range 0.3≤a2≤0.6 and fixing a0=2, a1=0.9.For case (I). The bifurcation diagram of (5) in  a0-xn and a0-yn space for a1=0.9, a3=0.01 is given in Figures 1(a) and 1(b) to show the dynamical changes of susceptible and infective, respectively, as  a0 varies. There is obvious phenomenon of bifurcation when we select the stepsize to be 10-4. The spectrum of Lyapunov exponents of the system (5) with respect to parameter a1 is given in Figure 1(c).(a) Bifurcation diagram in(a0,xn,yn) space. (b) The bifurcation diagram in a0-xn plane and in a0-yn plane. (c) Spectrum of Lyapunov exponents corresponding to (b). (d) The attractor of (5) with a0=2 (e) Time series of xn and yn of the same parameters.
(a)
(b)
(c)
(d)
(e)Moreover, we can see that the orbit with initial values(1,1) approaches to the stable fixed point P2 for a0<1.8 approximately, and Hopf bifurcation occurs at a0≈1.8. When increased to a0≈2.05, (5) becomes stable. In Figures 1(b) and 1(c), we observe the period -4,8 windows within the chaotic regions and boundary crisis at a0≈2.05. For a0∈(1.95,2.05)the maximum Lyapunov exponents are positive which correspond to chaotic region. To well see the dynamics, the attractor in the system (5) and time series of xn and yn are given in Figures 1(d) and 1(e), respectively.For case (II). Figure 2(a) is the bifurcation diagram of (5) for xn and yn and depicts that there are period -4,5 windows within the chaotic regions and boundary crisis at a1≈1. Figure 2(b) shows the spectrum of Lyapunov exponents of the system (5) with respect to parameter a1. For a1∈(0.85,0.95) the maximum Lyapunov exponents are positive which correspond to chaotic region.(a)
(b)For case (III). Figures 3(a) and 3(b) are the bifurcation diagrams of (5) for xn and yn, respectively, and depict that there are period -2,4 windows within the chaotic regions and boundary crisis at a2=0.5. Figure 3(c) shows the spectrum of Lyapunov exponents of (5) with respect to parameter a2. Figure 3(d) is the attractor of (5) with a2=0.46.(a)
(b)
(c)
(d)


## 4. Conclusion

In this paper, we investigate the behaviors of a discrete-time SIS epidemic model with nonlinear incidence rate, and find many complex and new interesting dynamical phenomena. Without the recovery rate of infectious hosts, (1) becomes the SI model (see [16]). Our theoretical analysis and numerical simulations have demonstrated that the model exhibits the variety of dynamical behaviors, which includes that the discrete epidemic model undergoes transcritical bifurcation, flip bifurcation, Hopf bifurcation, and chaos. Furthermore, chaos can cause the population to run a higher risk of extinction due to the unpredictability [24, 25]. Thus, how to control chaos in the epidemic model is very important, which needs further consideration.



---

*Source: 705601-2013-05-08.xml*