# Analysis of a Delayed SIR Model with  Nonlinear Incidence Rate

Authors: Jin-Zhu Zhang; Zhen Jin; Quan-Xing Liu; Zhi-Yu Zhang
Journal: Discrete Dynamics in Nature and Society
 (2008)

---

## Abstract

An SIR epidemic model with incubation time and
saturated incidence rate is formulated, where the susceptibles are
assumed to satisfy the logistic equation and the incidence term is
of saturated form with the susceptible. The threshold valueℜ0 determining whether the disease dies out is found. The
results obtained show that the global dynamics are completely
determined by the values of the threshold value ℜ0 and time
delay (i.e., incubation time length). If ℜ0 is less than one,
the disease-free equilibrium is globally asymptotically stable and
the disease always dies out, while if it exceeds one there will be
an endemic. By using the time delay as a bifurcation parameter, the
local stability for the endemic equilibrium is investigated, and the
conditions with respect to the system to be absolutely stable and
conditionally stable are derived. Numerical results demonstrate that
the system with time delay exhibits rich complex dynamics, such as
quasiperiodic and chaotic patterns.


---

## Body

## 1. Introduction

Epidemic models with
nonlinear incidence  have been studied by many authors, and related literature
of SIR disease transmission model is quite large, whereS denotes the number of
individuals that are susceptible to infection, I denotes the number of individuals that are
infectious,  and R denotes the number of individuals that
have been removed with immunity. For example, a detailed dynamical analysis of the nonlinear
incidence rate βIpSq (where β is the average number of contacts per
infective per unit time) is given by Liu et al. [1, 2], Hethcote and van den Driessche [3], Moghadas and Alexander [4],
Korobeinikov and Maini [5], and others. After a study of the cholera epidemic
spread in Bari in 1973, Capasso and Serio [6] introduced a saturated incidence rate g(I)S into epidemic models, the incidence rate seems
more reasonable than βIpSq because the number of effective contacts
between infective individuals and susceptible individuals may saturate at high
infective levels due to crowding of infective individuals or due to the
protection measures by the susceptible individuals [7], such as the incidence rate β(IpS/(1+αIq)) (see [1, 8] and references therein).Zhang and Chen [9] investigated a class of SIR epidemiological models
under assumption that the susceptible satisfies the logistic equation and the
incidence rate is of the form βISq.
More recently, Zhang et al. [10], serve as an extended version to [9], have carried out a
long-term qualitative analysis incorporating incubation time delay into incidence
rate in the case of q=1,
that is, with the force of infection βS(t)I(t−τ),
which was proposed by Cooke [11]. The incubation period τ(τ>0) is a time, during which the infectious agents
develop in the vector, and only after that time the infected vector becomes
itself infectious. The detailed biological meanings and transmission mechanisms
were given in [11].
The results obtained in [10]
represent that, if the epidemic is persistent, introducing time delay
changes the dynamical behaviors of the epidemic state. If the threshold value
determining whether the disease dies out is larger than one and less than
three, the endemic equilibrium is absolutely stable (in the sense that it is
asymptotically stable for all values of the delays [12]); when it exceeds three,
the endemic equilibrium is conditionally stable (i.e., it is asymptotically
stable for the delays in some intervals), and limit cycles arise by Hopf-type
bifurcation with increasing time delay.In 1978, May and Anderson [13] proposed the saturated
incidence rate of the form β(SI/(1+αS)),
and used by some authors [14–17], recently. The effect of saturation factor (refer to α) stems from epidemical control. In the
absence of effective therapeutic treatment and vaccine, the epidemical control
strategies are based on taking appropriate preventive measures. For example, if
transmission vector is mosquito, these measures include mosquito reduction
mechanisms and personal protection against exposure to mosquitos. Mosquito
reduction mechanisms entail the elimination of mosquito breeding sites (such as
clearing culverts, roadside ditches, eliminating standing water, etc.),
larvaciding (killing of larvae before they become adults) and adulticiding
(killing of adult mosquitoes by spraying). On the other hand, personal
protection is based on preventing vector mosquitoes from biting humans (by
using mosquito repellents, avoiding locations where mosquitoes are biting, and
using barrier methods such as window screens and long-sleeved clothing)
[18–20].From a practical point of view, instead of the
bilinear incidence rate in [10], we consider saturation incidence rate in this paper
and assume the force of infection is in this version β(SI(t−τ)/(1+αS)) which is saturated with the susceptible. The
susceptible host population is also assumed to have the logistic growth with
carrying capacity K, with a specific growth rate constant r.
We can get a generalized SIR epidemiological model as follows:(1.1)S˙(t)=r(1−S(t)K)S(t)−βS(t)1+αS(t)I(t−τ),I˙(t)=βS(t)1+αS(t)I(t−τ)−μ1I(t)−γI(t),R˙(t)=γI(t)−μ2R(t),where K,r,α,γ,μ1, and μ2 are positive constants. α is the parameter that measures the inhibitory
effect, γ is the natural recovery rate of the infective
individuals, μ1 and μ2 represent the per capita death rates of
infectives and recovered, respectively. Notice that when α=0,
the system (1.1) becomes the system of bilinear incidence rate in [10], throughout this paper, we
assume α≠0.By mathematical analysis, we derive a threshold valueℜ0 and prove that the values of ℜ0 and incubation time length completely
determine the global dynamics of system (1.1), that is, this two factors
determine whether the disease approaches an endemic value or whether solutions
oscillate. If ℜ0≤1,
the disease-free equilibrium is globally asymptotically stable and the disease
always dies out, whereas if ℜ0>1,
the disease persists if it is initially present. By taking the incubation time
delay as a bifurcation parameter, the local stability for the endemic
equilibrium is investigated, and the conditions with respect to the system to
be absolutely stable and conditionally stable are derived. Numerical
simulations show that the system with time delay admits rich complex dynamic,
and a sequence of periodic solutions will emanate with increasing time delay,
which exhibits quite complex periodic and chaotic patterns.We arrange our paper as follows. In Section2, results
on positivity and boundedness of solutions are presented. In addition, we also
consider the equilibria of system (1.1) and give the threshold for the
existence of endemic equilibrium. In Section 3, we consider the global
stability of the disease-free equilibrium and obtain the necessary and sufficient
conditions for the permanence of endemic equilibrium. The local stability
analysis of system (1.1) is considered in Section 4. Some numerical results
will be given as applications in Section 5.


## 2. Preliminary Results

The initial conditionsϕ=(ϕ1,ϕ2,ϕ3) of (1.1) are
defined in the Banach space(2.1)C+={ϕ∈C([−τ,0],R+3):ϕ1(θ)=S(θ),ϕ2(θ)=I(θ),ϕ3(θ)=R(θ)},where R+3={(S,I,R)∈R3:S≥0,I≥0,R≥0}.
By a biological meaning, we assume that ϕi(0)>0(i=1,2,3).It can be verified that the positive coneR+3 is positively invariant with respect to (1.1)
from [21, Lemma 2.1].Lemma 2.1.
All
feasible solutions of the system (1.1) are bounded and enter the region(2.2)Ωε={(S,I,R)∈R+3:S+I+R≤(r+1)μmM+ε,∀ε>0},where μm=min{1,μ1,μ2},limsupt→∞S(t)≤M:=max{S(0),K}.Proof.
From the first
equation of (1.1), we get(2.3)S˙(t)≤r(1−S(t)K)S(t),by comparison, we have limsupt→∞S(t)≤M. The total host population size N(t) can be determined by N(t)=S(t)+I(t)+R(t),
and(2.4)N˙(t)=r(1−S(t)K)S(t)−μ1I(t)−μ2R(t)≤(r+1)S(t)−S(t)−μ1I(t)−μ2R(t)≤(r+1)M−μmN.Thus, we have 0≤N≤((r+1)/μm)M, as t→∞.
Therefore, all feasible solutions of the system (1.1) are bounded and enter the
region Ωε.
This completes the proof of Lemma 2.1.Lemma2.1 shows that the solutions of system (1.1) are
bounded and, hence, lie in a compact set and are continuable for all positive
time.Letℜ0=K[β−α(μ1+γ)]/(μ1+γ).
For system (1.1), there always exists the equilibria E0=(0,0,0), E1=(K,0,0),
if ℜ0>1,
there also exists an endemic equilibrium E+=(S*,I*,R*),
where(2.5)S*=μ1+γβ−α(μ1+γ),I*=rS*2K(μ1+γ)(ℜ0−1),R*=γμ2I*.


## 3. Permanence

Before starting our theorem, we give the following
lemma.Lemma 3.1 (see[22]).
Consider
the following equation:(3.1)u˙(t)=au(t−τ)−bu(t),where a,b,τ>0;
and u(t)>0 for −τ≤t≤0.
One has
(i)
ifa<b,
then limt→∞u(t)=0;(ii)
ifa>b,
then limt→∞u(t)=+∞.Theorem 3.2.
Ifℜ0≤1,
then the solutions of (1.1), with respect to Ωε for any ε,
satisfy (S(t),I(t),R(t))→(K,0,0) as t→∞.Proof.
We consider
first the case whenℜ0<1.
From the first equation of (1.1), then there exists a ε>0 such that S(t)<K+ε for some T1>0 when t≥T1. Since S/(1+αS) is increasing function with respect to S,
then from the second equation of (1.1), we have(3.2)I˙(t)≤βK+ε1+α(K+ε)I(t−τ)−(μ1+γ)I(t).
Sinceℜ0<1, we have(3.3)βK+ε1+α(K+ε)<μ1+γ.By Lemma 3.1, we have limsupt→∞I(t)=0.limsupt→∞S(t)=K in terms of S(t)=K is the global attractivor of S(t)=r(1−S(t)/K)S(t).
Next, we shall consider the case whenℜ0=1. Noticing that ℜ0=1 is equal to β(K/(1+αK))=μ1+γ. Since S′(t)≤r(K−S(t))S(t), S(t) is always decreasing when above K.
If S(t) should ever get below K then S(t) must stay strictly below K for all subsequent time. This implies there
are two possible cases, either
(i)
S(t)→K from above as t→∞,
or(ii)
there existsT such that S(t)<K for all t>T.
In the first of these cases, we have only to show thatI(t)→0. Integrating the first equation for S from τ to t+τ in (1.1), we get(3.4)S(t+τ)−S(τ)=∫τt+τrS(u)(1−S(u)K)du−∫τt+τβS(u)1+αS(u)I(u−τ)du,≤∫τt+τrS(u)(1−S(u)K)du−∫τt+τβK1+αKI(u−τ)du,=∫τt+τrS(u)(1−S(u)K)du−∫τt+τ(μ1+γ)I(u−τ)du.Then,(3.5)(μ1+γ)∫0tI(u)du≤∫τt+τrS(u)(1−S(u)K)︸<0du−S(t+τ)+S(τ)≤S(τ)≤S(0).Letting t→∞,
we conclude that I(t)∈L1(0,∞) and, therefore, I(t)→0.
In the second of these cases, consider the
functional(3.6)V=I(t)+(μ1+γ)∫t−τtI(u)du.Then, for all t>T+τ,(3.7)V˙(t)|(1.1)=I˙(t)+(μ1+γ)[I(t)−I(t−τ)]=βS(t)1+αS(t)I(t−τ)−(μ1+γ)I(t−τ)=β[S(t)1+αS(t)−K1+αK]I(u−τ)<0.
A direct application of Liapunov-LaSalle type theorem
[22] shows that limt→∞I(t)=0.
By the third equation of (1.1), we get that limt→∞I(t)=0 implies limt→∞R(t)=0. This proves ℜ0≤1 is the sufficient condition for limt→∞(S(t),I(t),R(t))=(K,0,0).According to (1.1) and the definitions for permanence
in [23], we have the
following lemma.Lemma 3.3.
Permanence
ofS(t),I(t) in system (1.1) implies that of R(t).Next, we represent our main results in this section.Theorem 3.4.
System (1.1) is permanent if it satisfies ℜ0>1.In order to prove Theorem3.4, we present uniform
persistence theory for infinite dimensional systems from [24]. Let X be a complete metric
space. Suppose that X0 is open and dense in X and X0∪X0=X,X0∩X0=∅.
Assume that T(x) is a C0 semigroup on X satisfying(3.8)T(t):X0⟶X0,T(t):X0⟶X0.Let Tb(t)=T(t)|X0 and let Ab be the global attractor for Tb(t).Lemma 3.5 (see[24]).
Suppose
thatT(t) satisfies (3.8) and one has the following:
(i)
there
is at0≥0 such that T(t) is compact for t>t0;(ii)
T(t) is point dissipative in X;(iii)
Ab=⋃x∈Abω(x) is isolated and thus has an acyclic covering M^, where (3.9)M^={M1,M2,…,Mn};(iv)
Ws(Mi)∩X0=∅ for i=1,2,…,n.ThenX0 is a uniform repeller with respect to X0,
that is, there is an ε>0 such that for any x∈X0, liminft→+∞d(T(t)x,X0)≥ε,
where d is the distance of T(t)x from X0.Proof.
We first prove thatℜ0>1 leads to the permanence of system (1.1).
By Lemma3.3, we only need to consider the following
subsystem of (1.1) and prove that (S(t),I(t)) in system (3.10) are permanent if and only if ℜ0>1 holds(3.10)S˙(t)=r(1−S(t)K)S(t)−βS(t)1+αS(t)I(t−τ),I˙(t)=βS(t)1+αS(t)I(t−τ)−μ1I(t)−γI(t),where S(θ),I(θ)≥0 are continuous on −τ≤θ≤0,
and S(0),I(0)>0.
We begin by showing that the boundary planes ofℝ+2={(x,y):x≥0,y≥0} repel the positive solutions to system (3.10)
uniformly. Let C+([−τ,0],ℝ+2) denote the space of continuous functions
mapping [−τ,0] into ℝ+2.
We choose(3.11)C1={(φ1,φ2)∈C+([−τ,0],ℝ+2):φ1(θ)≡0,θ∈[−τ,0]},C2={(φ1,φ2)∈C+([−τ,0],ℝ+2):φ1(θ)>0,φ2(θ)≡0,θ∈[−τ,0]}.
Denote C0=C1∪C2,X=C+([−τ,0],ℝ+2) and C0= Int C+([−τ,0],ℝ+2).
Next, we verify that the conditions of Lemma3.5 are
satisfied. By the definition of C0 and C0 and system (3.10), it is easy to see that C0 and C0 are positively invariant. Moreover, conditions
(i) and (ii) of Lemma 3.5 are clearly satisfied. Thus, we only need to verify
conditions (iii) and (iv). Since system (3.10)
possesses two constant solutions in C0:E0∈C1,E1∈C2 with(3.12)E0={(φ1,φ2)∈C+([−τ,0],ℝ+2):φ1(θ)≡φ2(θ)≡0,θ∈[−τ,0]},E1={(φ1,φ2)∈C+([−τ,0],ℝ+2):φ1(θ)≡K,φ2(θ)≡0,θ∈[−τ,0]},and we have S˙(t)|(φ1,φ2)∈C1≡0,
then we get S(t)|(φ1,φ2)∈C1≡0 for all t≥0,
according to the second equation of (3.10), we have I˙(t)|(φ1,φ2)∈C1=−(μ1+γ)I(t)≤0,
hence all points in C1 approach E0,
C1=Ws(E0).
Similarly, we have all points in C2 approach E1,
that is, C2=Ws(E1).
This shows that invariant sets E0 and E1 are isolated invariant, then {E0,E1} is isolated and is an acyclic covering,
satisfying condition (iii) of Lemma 3.5.
Now, we show thatWs(Ei)∩C0=∅,i=0,1. We only need to prove Ws(E1)∩C0=∅,
since the proof for Ws(E0)∩C0=∅ is simple.
Assume the contrary, that is,Ws(E1)∩C0≠∅,
then there exists a positive solution (S(t),I(t)) to system (3.10) with limt→∞(S(t),I(t))=(K,0). Since ℜ0>1,
then for a sufficiently small ε>0 with μ1+γ<β((K−ε)/(1+α(K−ε))),
there exists a positive constant T=T(ε) such that (3.13)S(t)>K−ε>0,0<I(t)<ε∀t≥T.
By the second
equation of (3.10) and noting that β(S/(1+αS)) is an increasing function with respect to S,
then we have(3.14)I˙(t)≥βK−ε1+α(K−ε)I(t−τ)−(μ1+γ)I(t),t≥T+τ.According to the comparison
principle, limt→∞I(t)=∞ when ℜ0>1,
contradicting I(t)<ε.
Then we have Ws(E1)∩C0=∅.
At this time, we are able to conclude from Lemma 3.5 that C0 repels the positive solutions of (3.10)
uniformly. Incorporating the above results into Lemmas 3.3 and 3.5, we know
that system (1.1) is permanent.
Next, we verify that permanence of system (1.1)
indicates ℜ0>1.
Assume that the contrary holds, that is, ℜ0≤1,
then by Theorem 3.2, (S(t),I(t),R(t))→(K,0,0), as t→∞, contradicting permanence of system (1.1). This
proves Theorem 3.4.


## 4. Linearized Analysis

The dynamics of model (1.1) are determined by the
first two equations. Therefore, throughout the remainder of this paper, we
consider the subsystem (3.10), and rewrite it as follows:(4.1)S˙(t)=r(1−S(t)K)S(t)−βS(t)1+αS(t)I(t−τ),I˙(t)=βS(t)1+αS(t)I(t−τ)−μ1I(t)−γI(t).LetE^=(S^,I^) be any equilibrium of (4.1), linearized system
of (4.1) at E^=(S^,I^),
we get(4.2)x˙(t)=[r−2rKS^−βI^(1+αS^)2]x(t)−βS^1+αS^y(t−τ),y˙(t)=βI^(1+αS^)2x(t)+βS^1+αS^y(t−τ)−(μ1+γ)y(t).Then the characteristic equation
of (4.1) at E^ is given by(4.3)det[r−2rKS^−βI^(1+αS^)2−λ−βS^1+αS^e−λτβI^(1+αS^)2βS^1+αS^e−λτ−(μ1+γ)−λ]=0.At the equilibrium E^0=(0,0),
characteristic equation (4.3) reduces to(4.4)(λ−r)(λ+μ1+γ)=0.
Obviously,
(4.4) has a positive root λ=r independent of any parameters. Hence, E^0 is always a unstable saddle point.Theorem 4.1.
For the
system (4.1), the equilibrium E^1=(K,0) is(i)
asymptotic stable ifℜ0<1;(ii)
linearly
neutrally stable ifℜ0=1;(iii)
unstable
ifℜ0>1.Proof.
The characteristic equation atE^1 is(4.5)(λ+r)[λ+(μ1+γ)(1−ℜ0+αK1+αKe−λτ)]=0.
Equation (4.5) has a negative real part characteristic
root λ=−r and roots of(4.6)F(λ)=λ+(μ1+γ)(1−ℜ0+αK1+αKe−λτ)=0.
(i) Assume thatℜ0<1,
(4.6) has characteristic root λ=(μ1+γ)((ℜ0+αK)/(1+αK)−1)<0 when τ=0.
If λ=iω is a root of (4.6), it must
satisfy(4.7)ω2=(μ1+γ)2[(ℜ0+αK1+αK)2−1].When ℜ0<1,
there are no positive real roots ω.
This shows that all roots of F(λ)=0 must have negative real parts, therefore, E^1 is an asymptotically stable equilibrium.
(ii) Assume thatℜ0=1,
then λ=0 is a root of (4.6). It is easy to verify that λ=0 is a simple characteristic root. If the other
roots are λ=α+iω,
then they must satisfy(4.8)[α+(μ1+γ)]2+ω2=(μ1+γ)2e−2ατ,and we must have α≤0.
Therefore E^1 is linearly neutrally stable.
(iii) Assume thatℜ0>1,
then F(0)<0,
and F(+∞)=+∞.
Hence, F(λ) has at least one positive root and E^1 is unstable.By the arguments to Theorems3.2 and 4.1, we
directly have the following corollary.Corollary 4.2.
The equilibriumE1=(K,0,0) of system (1.1) is global asymptotically
stable if ℜ0≤1 holds true in the feasible region Ωε for any ε>0.In the following, we will study the linear stability
of the positive equilibriumE^+=(S*,I*) of (4.1). We can see that the characteristic
roots of (4.3) at positive equilibrium E^+ are the roots of(4.9)det[r(1−2ℜ0)−r1+αS^(1−1ℜ0)−λ−(μ1+γ)e−λτr1+αS^(1−1ℜ0)(μ1+γ)e−λτ−(μ1+γ)−λ]=0.
Since β(I*/(1+αS*))=r(1−1/ℜ0) and β(S*/(1+αS*))=μ1+γ at (S*,I*),
we have(4.10)P(λ,τ)+Q(λ,τ)e−λτ=0,where(4.11)P(λ,τ)=λ2+λ[−r(1−2ℜ0)+(μ1+γ)+r1+αS*(1−1ℜ0)]+(μ1+γ)[−r(1−2ℜ0)+r1+αS*(1−1ℜ0)],Q(λ,τ)=−λ(μ1+γ)+r(μ1+γ)(1−2ℜ0).
When τ=0,
the DDE (4.1) becomes ODE which has the same equilibria E^ as follows:(4.12)S˙(t)=r(1−S(t)K)S(t)−βS(t)1+αS(t)I(t),I˙(t)=βS(t)1+αS(t)I(t)−μ1I(t)−γI(t),and (4.10)
becomes(4.13)λ2+λ[−r(1−2ℜ0)+r1+αS*(1−1ℜ0)]+(μ1+γ)r1+αS*(1−1ℜ0)=0.Defineℜcc=2+1/αS*.If−r(1−2/ℜ0)+(r/(1+αS*))(1−1/ℜ0)>0,
that is, 1<ℜ0<ℜcc=2+1/αS*,
the system (4.12) is locally asymptotically stable. If ℜ0>ℜcc,
the unique positive equilibrium of (4.12) is unstable, system (4.12) becomes
oscillatory in a stable limit cycle [25], and this limit cycle is unique [26, 27].Ifλ=iω(ω>0) is a root of (4.10), then by separating the
real and imaginary parts, we get(4.14)−ω2+(μ1+γ)[−r(1−2ℜ0)+r1+αS*(1−1ℜ0)]=−r(1−2ℜ0)(μ1+γ)cosωτ+ω(μ1+γ)sinωτ,ω[−r(1−2ℜ0)+(μ1+γ)+r1+αS*(1−1ℜ0)]=ω(μ1+γ)cosωτ+r(1−2ℜ0)(μ1+γ)sinωτ.
Squaring and
adding both equations, then we have (4.15)ω4+ω2[−r(1−2ℜ0)+r1+αS*(1−1ℜ0)]2+r1+αS*(1−1ℜ0)(μ1+γ)2[−2r(1−2ℜ0)+r1+αS*(1−1ℜ0)]=0.Defineℜc=2+1/(1+2αS*).If−2r(1−2/ℜ0)+(r/(1+αS*))(1−1/ℜ0)≥0,
that is, 1<ℜ0≤ℜc,
there is no positive real ω satisfying (4.15), thus eigenvalues of (4.10)
do not approach the imaginary axis for any τ>0.
This shows that E^+ is absolutely stable when 1<ℜ0≤ℜc.Ifℜ0>ℜc,
there is a unique positive ω0 satisfying (4.15). That is, (4.10) has a unique
pair of purely imaginary roots ±iω0.From (4.14), τn corresponding to ω0 can be obtained as follows:(4.16)τn=1ω0arccos{ω02ℤ−ℂ(μ1+γ)r(1−2/ℜ0)ω02(μ1+γ)+[r(1−2/ℜ0)]2(μ1+γ)}+2nπω0,n=0,1,2,…,
where ℤ denotes [(r/(1+αS*))(1−1/ℜ0)+(μ1+γ)] and ℂ denotes [−r(1−2/ℜ0)+(r/(1+αS*))(1−1/ℜ0)].Further,(4.17)dRe(λ)dτ|λ=iω0=Re(dλdτ)−1|λ=iω0=Re{λ2−(μ1+γ)[−r(1−2/ℜ0)+(r/(1+αS*))(1−1/ℜ0)]−λ2P(λ,τ)}λ=iω0+Re[−r(1−2/ℜ0)(μ1+γ)λ2Q(λ,τ)]λ=iω0=ω04−(μ1+γ)2[−r(1−2/ℜ0)+(r/(1+αS*))(1−1/ℜ0)]2+[r(1−2/ℜ0)(μ1+γ)]2(μ1+γ)2+[r(1−2/ℜ0)(μ1+γ)]2=ω04−(μ1+γ)2(r/(1+αS*))(1−1/ℜ0)[−2r(1−2/ℜ0)+(r/(1+αS*))(1−1/ℜ0)](μ1+γ)2+[r(1−2/ℜ0)(μ1+γ)]2.
Under the
condition ℜ0>ℜc,
that is, −2r(1−2/ℜ0)+(r/(1+αS*))(1−1/ℜ0)<0,
then we have
dRe(λ(τ))/dτ∣λ=iω0>0.Ifℜc<ℜ0<ℜcc,
there exists a critical value τ0, when τ<τ0,E^+ is stable; when τ>τ0,E^+ is unstable.Summarizing the discussion above, we have the
following conclusion.Theorem 4.3.
For system (4.1), one has(i)
if1<ℜ0≤ℜc holds true, then E^+ is absolutely stable;(ii)
ifℜc<ℜ0<ℜcc holds true, then E^+ is conditionally stable, that is, there is a
critical delay value τ0 such that E^+ is asymptotically stable when τ∈[0,τ0) and unstable when τ>τ0.
Furthermore, system (4.1) undergoes Hopf bifurcation at E^+ when τ=τn,n=0,1,2,…;(iii)
ifℜ0>ℜcc,
then there is also a critical delay value τ0 such that the periodic solution is still
stable when τ∈[0,τ0),
however, there are a sequence of periodic solutions emanate when τ=τn,n=0,1,2,….These results are summarized in Table1.Table 1
Compare DDE (4.1) with ODE (4.12).
Case0<ℜ0≤11<ℜ0≤ℜcℜc<ℜ0<ℜccℜ0>ℜccODEE^1 GASE^+ LASE^+ LASstable periodic solutionDDEE^1 GASE^+ ALASaE^+ CLASbcomplex dynamic phenomenaaabsolutely stable
bconditionally stableRemark 4.4.
In fact, in the case ofℜ0>ℜcc,
we have known that ODE (4.12) has a stable periodic solution [25]. If we consider the impact
of the incubation time on ODE (4.12), that is, DDE (4.1), from above discussion,
we can see that there are a sequence of periodic solutions bifurcate from the
positive equilibrium E^+ when the time delay takes the critical delay τn such that previous stable periodic solution
losses stability, which will lead to complex dynamic phenomena. This can be
seen from the simulation results Section 5.In addition, we want to mention that Theorem4.3(ii)
cannot determine the direction and stability of bifurcation periodic solutions,
this can be done by analyzing the high-order terms in terms of [28]. The method is very complex
and trivial, here we omit it.From the point of biology, in comparison with the
results of [10], we
see that if there is the inhibition effect from the behavioral change of the
susceptible individuals when the infective increases (i.e., we take use of the
saturation incidence rate), then the threshold ℜc decline and less than 3. It is crucial for the
government to take the corresponding control measures and policies against the
disease when the epidemic outbreaks.


## 5. Numerical Results

The main goal of the previous section was to
qualitatively characterize the dynamic behaviors of system (1.1) at long term.
In this section, we confirm our previous theoretical analysis in Section 4 and
demonstrate that the local behaviors in the regions of the parameters space
correspond to complex population dynamics to system (1.1). The objective is to
explore the possibility of chaotic behavior in system (1.1). It is difficult to
test whether there exists chaos in a time-delayed system, but numerical
simulation analysis is a valid method for such a system. Extensive numerical
simulations are carried out for different values of saturation parameter α and recover rate γ.
The quality results are as follows.First, we consider the property of system (1.1) in the
regions of the parameter space corresponds to complex population dynamics in
the case of ℜ0>ℜcc.
To illustrate the transition from the periodic pattern to chaotic pattern, we
concentrate on the regions of small and large time delay as an example. We
consider the set of parameter values as r=0.1, K=80, β=0.1, μ1=0.5, μ2=0.1, α=0.05, and γ=0.5.
By calculating, there exists the relation of ℜ0>ℜcc (ℜ0=4.00, ℜcc=3.00) for system (1.1). Then, from Theorem
4.3(iii), the system (1.1) has a stable period solution if the time delay is
less than the critical delay τ0≐0.23 (see Figure 1), and the periodic solution will
lost stability when the time delay is greater than the critical delay τ0≐0.23,
and then a typical chaos was observed with increasing the time delay (see
Figure 2). This phenomenon has been verified by the bifurcation diagram via
delay τ,
as shown in Figure 3.Temporal behavior of the
infected and corresponding three-dimensional phase are plotted for the system
(1.1) subject to ℜ0>ℜcc (ℜ0=4.00, ℜcc=3.00). The parameters are r=0.1, K=80, β=0.1, μ1=0.5, μ2=0.1, α=0.05, γ=0.5, and τ=0.001 with initial value (28, 3, 7).
(a)(b)The parameters like in Figure1 but τ=30.
(a)(b)Figure 3
Bifurcation diagram of system (1.1):
successive maxima of the infected are plotted for increasing values of the time
delay τ,
with parameters r=0.1, K=80, β=0.1, μ1=0.5, μ2=0.1, α=0.05, and γ=0.5.We increase the recovery rateγ,
let γ=0.75,
and fix the other parameters as above. By calculating, the system (1.1)
satisfies ℜc<ℜ0<ℜcc (ℜc=2.23, ℜ0=2.40, ℜcc=2.6). In this context, by Theorem 4.3(ii), system
(1.1) is conditionally stable at unique positive equilibrium E+=(33.33,1.56,11.67),
that is, there exists a critical delay τ0≐0.31 such that E+ is stable if τ<τ0 (see Figure 4), E+ will lose stability by an Hopf bifurcation if τ>τ0,
as shown in Figure 5. We find that the periodic solution, quasiperiod, and
chaos patterns emerge with increasing time delay. This may be clear from the
bifurcation diagram (see Figure 7).Temporal behavior of
the infected and corresponding three-dimensional phase are plotted for the
system (1.1) subject to 1<ℜc<ℜ0<ℜcc (ℜc=2.23, ℜ0=2.40, ℜcc=2.6). The parameters are r=0.1, K=80, β=0.1, μ1=0.5, μ2=0.1, α=0.05, γ=0.75, and τ=0.001 with initial value (35, 10, 15).
(a)(b)Temporal behavior of the
infected and corresponding three-dimensional phase are plotted for the system
(1.1) subject to 1<ℜc<ℜ0<ℜcc (ℜc=2.23, ℜ0=2.40, ℜcc=2.6). The parameters are r=0.1, K=80, β=0.1, μ1=0.5, μ2=0.1, α=0.05, γ=0.75,
and τ=0.5 with initial value (40, 3, 7).
(a)(b)The parameters like in
Figure5 but τ=27.
(a)(b)Figure 7
Bifurcation diagram of system (1.1): successive maxima of the infected
are plotted for increasing time delay τ,
with parameters r=0.1, K=80, β=0.1, μ1=0.5, μ2=0.1, α=0.05m, and γ=0.75.Now, we fix the parameterγ=0.75,
change α, and let α=0.055 (i.e., we take some measures to protect on
susceptibles), and the other parameters are also as above. By calculating,
system (1.1) always exists a relationship 1<ℜ0<ℜc (ℜ0=2.00, ℜc=2.19). By Theorem 4.3(i), we know that system
(1.1) is absolutely stable at unique positive equilibrium E+=(40,1.60,12) for any value of the time delay, as shown in
Figure 8. Our numerical simulations have demonstrated the validity of our
theoretical analysis, that is, the values of threshold value ℜ0 and incubation time τ length completely determine the dynamics of
system (1.1).Temporal
behavior of the infected and corresponding three-dimensional phase are plotted
for the system (1.1) subject to 1<ℜ0<ℜc.
The parameters like in Figure 5 but α=0.055, γ=0.75,
and τ=1.5 with initial value (25, 4, 20).
(a)(b)It is necessary to indicate that the system (1.1) is
realistic at the initial phase of the disease emergence since the number of the
infected is rare. However, when the number of the infected is large, it is more
reasonable that one should replace the term r(1−S/K)S in the first equation of (1.1) by r(1−(S+I+R)/K)S.
Hence, a profound understanding for this case is still desirable and could
motivate further investigations.



---

*Source: 636153-2009-01-27.xml*