# Rich Dynamics of an Epidemic Model with Saturation Recovery

Authors: Hui Wan; Jing-an Cui
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

A SIR epidemic model is proposed to understand the impact of limited medical resource on infectious disease transmission. The basic reproduction number is identified. Existence and stability of equilibria are obtained under different conditions. Bifurcations, including backward bifurcation and Hopf bifurcation, are analyzed. Our results suggest that the model considering the impact of limited medical resource may exhibit vital dynamics, such as bistability and periodicity when the basic reproduction numberℝ0 is less than unity, which implies that the basic reproductive number itself is not enough to describe whether the disease will prevail or not and a subthreshold number is needed. It is also shown that a sufficient number of sickbeds and other medical resources are very important for disease control and eradication. Considering the costs, we provide a method to estimate a suitable treatment capacity for a disease in a region.


---

## Body

## 1. Introduction

In recent years, efforts have been made to develop realistic mathematical models for the transmission dynamics of known and emerging infectious diseases [1–8]. The development of such models aims at understanding the epidemiological transmission patterns and predicting the consequences of the introduction of public health interventions to control the possible outbreak and spread of the disease.In this paper, we will focus on the impact of limited medical resource or treatment capacity on infectious disease transmission. In fact, each city should have its suitable treatment capacity. If it is too large, the city pays for unnecessary costs. While if it is too small, the city has the risk of the outbreak of a disease. Statistics [9] show that, for example, there are 16.4 sickbeds for every 1000 residents in Japan on average in 1999, but the corresponding data are 2.6, 2.4, and 1.1, respectively, in Turkey, China, and Mexico. Thus, it is important to understand the impact of limited medical resource on disease transmission and to determine a suitable treatment capacity for a disease.In classical disease transmission models, the recovery from infected class per unit of time is assumed to be proportional to the number of infective individuals. However, it is not a reasonable approximation to the truth when the number of the infectious individuals is large and the treatment capacity of hospitals is researched, provided that the infected individuals cannot recover unless they were given timely treatment in hospitals. Then the number of patients needing to be treated in hospitals may exceed the number of the hospital beds. In this case, the recovery rate from infective class will be saturated at a maximum, especially in the rural areas in many developing countries.In paper [10], Wang and Ruan adopt a constant treatment rate, which simulates a limited capacity for treatment. Note that a constant treatment is suitable when the number of infectives is large; hence in [7], Wang modified the rate of treatment to another function which is proportional to the number of the infectives when the capacity of treatment is not reached and, otherwise, takes the maximal capacity. In order to model the impact of limited medical resource on infectious disease transmission precisely and provide a method to estimate a suitable treatment capacity for a disease, we propose an epidemic model with saturation recovery.The organization of this paper is as follows. In Section2, we propose an SIR epidemic model and analyze the disease-free equilibrium. In Section 3, we mainly discuss the existence of endemic equilibrium (EE). In Section 4, the stability of EE and bifurcation are presented. Section 5 is the discussion.


## 2. The Model

In order to give a more realistic recovery rate for epidemic models and study the impact of limited medical resource on disease transmission, in this paper, we adopt the Verhulst-type function ([11])
(1)h(I)=cIb+I
to model the treated part which is increasing for small infectives and approaches the maximum for large infectives. Hence c, as the limit of h(I) as I tends to infinity, is the maximum of the treatment capacity in a region, and b, the infected size at which is 50% saturation (h(b)=c/2), measures how soon saturation occurs. Here h(I) satisfies h(0)=0,h(∞)=c and dh/dI>0, which means that the treatment rate is increasing for a small number of infectives and approaches a maximum for large number of infectives.We classify the population in a given region/area into three categories: susceptible, infective, and recovered. LetS(t),I(t), and R(t) denote the number of susceptible, infective, and recovered individuals at time t, respectively. Provided that the infected individuals cannot recover unless they were given timely treatment in hospitals, based on standard SIR model with mass action incidence, we can construct a model
(2)dSdt=A-dS-βSI,dIdt=βSI-(d+ν)I-cIb+I,dRdt=cIb+I-dR,
where all the parameters are positive and(i)
Ais the recruitment rate of susceptible population;(ii)
d is natural death rate and ν is the disease-induced death rate;(iii)
c is the maximum of treatment per unit of time, and b, the infected size at which is 50% saturation (h(b)=c/2), measures how soon saturation occurs;(iv)
β is the transmission rate.Note that the first two equations are independent of the third one; we need only to study the following reduced model:(3)dSdt=A-dS-βSI,dIdt=βSI-(d+ν)I-cIb+I.Model (3) has one disease-free equilibrium at E0=(A/d,0). Using the formulae in [12], a straightforward calculation gives the reproduction number:
(4)ℝ0=Abβd(bd+bν+c).The disease-free equilibriumE0 has two eigenvalues -d and ℝ0-1. Therefore we have the following proposition.Proposition 1.
For the model (3), the disease-free equilibrium E0 is locally asymptotically stable if ℝ0<1 and unstable if ℝ0>1.


## 3. Existence of the Endemic Equilibrium

Let(5)h(I)=Ad+βI,g(I)=d+νβ+cβ(b+I).
Then we can rewrite the model (3) as
(6)dSdt=-(d+βI)(S-h(I)),dIdt=βI(S-g(I)).Let the right hand side of (6) be zero. If an endemic equilibrium exists, its (S,I) coordinates must satisfy
(7)S=g(I),S=h(I).
We note that limI→∞h(I)=0,limI→∞g(I)=(d+ν)/β,dh/dI<0,dg/dI<0, and d2h/dI2>0,d2g/dI2>0. In addition, ℝ0=h(0)/g(0).EliminatingS from (7) gives an equation of the form
(8)I2+a1I+a2=0,
where
(9)a1=(bβ+d)(d+ν)+cβ-Aββ(d+ν),a2=bd(d+ν)+cd-βAbβ(d+ν)=bd(d+ν)+cdβ(d+ν)(1-ℝ0).Ifℝ0>1, then a2<0, and there is a unique positive root for (8) which implies that a unique endemic equilibrium E*(S*,I*) exists (see Figure 1(b)).Existence and number of endemic equilibria.
(a)
No positive equilibrium
(b)
Unique positive equilibrium
(c)
Unique equilibrium of multiplicity 2
(d)
Two positive equilibriaIfℝ0=1, then a2=0 and there is a unique nonzero solution of (8) I=-a1 which is positive if and only if a1<0. Then, there is a unique endemic equilibrium E*(S*,I*) when a1<0, and there are not endemic equilibria when a1≥0.Ifℝ0<1, then a2>0. For (8) to have at least one positive root we must have
(10)a1<0,Δ≜a12-4a2⩾0.
Solving Δ=0 in terms of ℝ0, one gets ℝ0=ℝ^0 where
(11)ℝ^0(c)=Ab(d+ν)[b(d+ν)+(A+c)2]×((A-c)2(c+b(d+ν))[(A-c)2b(d+ν)(b(d+ν)+2(A+c))+(A-c)2])-1.
One can verify that, provided a1<0, model (3) has exactly 0, 1, and 2 endemic equilibria as shown in Figures 1(a), 1(c), and 1(d) for ℝ0<ℝ^0,ℝ0=ℝ^0,andℝ0>ℝ^0, respectively.In summary, regarding the existence of endemic equilibrium, we have the following.Theorem 2.
For model (3), one has the following.(1)
Ifℝ0>1, there exists a unique positive equilibrium E*(S*,I*).(2)
Ifℝ0=1; there is a positive equilibrium E*(S*,I*) when a1<0; otherwise there is no positive equilibrium.(3)
Ifℝ0<1 and a1≥0, there is no positive equilibrium.(4)
Ifℝ^0<ℝ0<1 and a1<0, there are two positive equilibria E* and E*.(5)
Ifℝ0=ℝ^0 and a1<0,E* and E* coalesce together as a unique equilibrium of multiplicity two.(6)
Ifℝ0<ℝ^0 and a1<0, there is no positive equilibrium.
When exist,E*(S*,I*) and E*(S*,I*) are the corresponding equilibria, and
(12)S*=h(I*),S*=h(I*),I*=-a1+Δ2,I*=-a1-Δ2.c, the maximum treatment per unit of time, is related to the treatment capacity in a region and the recruitment rate of susceptible individuals A can be changed by vaccination, and so forth; they are all important for disease control. In the following section, we choose c and A as parameters to discuss the distribution of the equilibria of system (6) in plane (c,A). Solving equations a1=0 and a2=0 for A, we can get A=A1(c) and A=A2(c), respectively, where
(13)A1(c)=c+(b+dβ)(d+ν),A2(c)=dbβc+d(d+ν)β.
Then
(14)a1=(A1(c)-A)1d+ν,a2=(A2(c)-A)bd+ν.Obviously,a1>0 (or a1<0) when A<A1(c) (or A>A1(c)) and a2>0 (or a2<0) when A<A2(c) (or A>A2(c)).It is easy to see that ifA>A2(c), (7) has a unique positive solution.Now we consider the case whenA⩽A2(c) (the case a2≥0). Writing the discriminant Δ=a12-4a2 as a function of c and A, we have
(15)Δ(c,A)=[(bβ+d)(d+ν)+cβ-Aββ(d+ν)]2-4[bd(d+ν)+cd-βAbβ(d+ν)].
If Δ=0 defines curves in the first quadrant of the (c,A) plane, then it must satisfy A⩽A2(c). Furthermore, if (8) has nonnegative roots, we need A1(c)<A⩽A2(c) and Δ⩾0. Calculate the difference of A2(c) and A1(c):
(16)A2(c)-A1(c)=d-bβbβc-b(d+ν).Ifd-bβ⩽0, then A2(c)<A1(c) and system (6) does not have any endemic equilibrium when A⩽A2(c).Ifd-bβ>0, one easily gets that c=c* from A2(c)-A1(c)=0, where
(17)c*=b2β(d+ν)d-bβ.
We have
(18)c⩽c*⟺A1(c)⩾A2(c),c>c*⟺A1(c)<A2(c).
Therefore, when c⩽c*, system (6) does not have any endemic equilibrium in this case. If c>c*, then system (6) has two equilibria when A1(c)<A<A2(c) and Δ>0; system (6) has a unique endemic equilibrium when A=A2(c); system (6) has a unique endemic equilibrium of multiplicity 2 when A1(c)<A<A2(c) and Δ=0; and system (6) has no equilibrium when Δ<0.Solving the equationΔ=0 for A, one gets two curves
(19)Γ1:A=A3(c),Γ2:A=A4(c),
in the (c,A) plane, where
(20)A3(c)=(d-bβ)(d+ν)+cβ+2cβ(d-bβ)(d+ν)β,A4(c)=(d-bβ)(d+v)+cβ-2cβ(d-bβ)(d+ν)β.
For any point of (c,A)∈Γ2, we have A<A1(c) by simple calculation. The two lines A=A1(c),A=A2(c) joint the curve Γ1:A=A3(c) at a point P(A1(c*),c*) when c=c*. For any c>c*, the curve Γ1 lies above A=A1(c) and under A=A2(c) as shown in Figure 2. In fact, we have Δ<0 for the point (c,A) on A=A1(c) and Δ>0 for the point (c,A) on A=A2(c), and hence the curve Γ1 lies between the two lines A=A1(c) and A=A2(c) when c>c*. As is shown in Figure 2, in the case when d-bβ>0, the curve A=A2(c) and the curve segment Γ1 subdivide the positive cone of the parameters plane (c,A) into three regions D0,D1 and D2, where
(21)D0={(c,A)∣0<c⩽c*,0<A⩽A2(c)}∪{(c,A)∣c>c*,0<A<A3(c)},D1={(c,A)∣0<c⩽c*,A>A2(c)}∪{(c,A)∣c>c*,A⩾A2(c)},D2={(c,A)∣c>c*,A3(c)<A<A2(c)}.
For parameter (c,A) in each region, the system has exactly 0, 1, and 2 positive equilibria, respectively. In summary, we have the following theorem regarding the number of endemic equilibria.Figure 2
The distribution of equilibrium on the plane of(c,A) when d-bβ>0. There are 0, 1, 2 equilibria in D0,D1,D2, respectively.Theorem 3.
For the model (6), with c*,A1(c),A2(c), and A3(c) defined as previously mentioned,(1)
ifd-bβ>0, as shown in Figure 2, one has the following.
(a)
For(c,A)∈D0, the model (6) has no equilibrium.(b)
For(c,A)∈D1, the model (6) has a unique equilibrium E*.(c)
For(c,A)∈D2, the model (6) has two equilibria E* and E*.(d)
For(c,A)∈Γ1, when c>c*,E* and E* coalesce at a unique endemic equilibrium of multiplicity 2.(e)
At the pointP(A1(c*),c*), the model (6) has no endemic equilibrium; the disease-free equilibrium is of multiplicity 2.(2)
On the other hand, ifd-bβ⩽0, then
(a)
ifA>A2(c), there is a unique endemic equilibrium E*;(b)
ifA⩽A2(c), there is no endemic equilibrium.
When exist,E*(h(I*),I*) and E*(h(I*),I*) are the corresponding equilibria, and
(22)I*=12[-a1+Δ(c,A)],I*=12[-a1-Δ(c,A)].Now we can prove the global stability of the disease-free equilibriumE0 when (c,A)∈D0 according to Theorem 3.Corollary 4.
If(c,A)∈D0, then the disease-free equilibrium E0 of (6) is globally asymptotically stable.Proof.
LetΩ={(S,I)∣S,I≥0,S+I≤A/d}. We will prove that Ω is a positively invariant set of (6). By (6), we have (dS/dt)|S=0=A>0 and (dI/dt)|I=0=0. Denote N(t)=S(t)+I(t); then
(23)dNdt|N=A/d≤A-dN|N=A/d=0.
Hence Ω is a positively invariant set of (6), and it attracts all the positive orbits of (6) state in R2+. We have shown that (6) does not have any positive equilibria when (c,A)∈D0. It follows from the Poincaré-Bendixson theorem that no periodic orbits (limit cycle) exist in Ω. Since Ω is a bounded positively invariant region for the model and E0 is the unique equilibrium in Ω, the local stability of E0 implies that the ω-limit set of every solution with an initial point in R2+ is E0. Hence, the disease-free equilibrium E0 is globally asymptotically stable.


## 4. Stability and Bifurcations of the Endemic Equilibria

In this section, we study the stability and bifurcation of the endemic equilibrium. Evaluating the Jacobian of the model (6) at E*(S*,I*) gives
(24)J=(-d-βI*-βS*βI*cI*(b+I*)2).
Then the characteristic equation about E* is given by
(25)λ2+H(I*,c)λ+I*G(I*,c)=0,
where
(26)H(I*,c)=d+βI*-cI*(b+I*)2,G(I,c)=Aβ2d+βI*-c(d+βI*)(b+I*)2.


### 4.1. Case 1:ℝ0>1

Whenℝ0>1, the model (6) has a unique endemic equilibrium E*. From (24) to (26), we have the following.Theorem 5.
Supposeℝ0>1. Then the disease endemic equilibrium E* of (6) is a stable node or focus when H(I*,c)>0;E* is an unstable node or focus when H(I*,c)<0, and (6) has at least one closed orbit in Ω;E* is a center of the linear system when H(I*,c)=0.Proof.
RewriteG(I,c) as
(27)G(I,c)=β(d+βI)[-cβ(b+I)2+βA(d+βI)2]=β(d+βI)[dg(I)dI-dh(I)dI].
Note that dg(I*)/dI-dh(I*)/dI>0; hence G(I*,c)>0. The unique endemic equilibrium E* is a stable node or focus when H(I*,c)>0,E* is an unstable node or focus when H(I*,c)<0;E* is a center of the linear system when H(I*,c)=0. If H(I*,c)<0, then E is an unstable node or focus. System (6) has at least one closed orbit in Ω from the Poincaré-Bendixson Theorem.As an example, we fixA=15000,β=0.0005,d=0.02,b=1,ν=34. When c=30, then ℝ0=5.86; the model (6) has unique endemic equilibrium E(S*,I*)=(68189.645,399.949) and H(I*,c)=0.14534>0. From Figure 3(a), we can see that the endemic equilibrium E* is stable. When c=126.99, then ℝ0=2.33; the model also has a unique endemic equilibrium E(S*,I*)=(68678.432,396.818) but H(I*,c)=-0.01<0. From Figure 3(b), we can see that the endemic equilibrium E* is unstable and there must be a stable limit cycle around E(S*,I*).The phase portraits of the model where the directions of trajectories are counterclockwise.
(a)
H
(
I
*
,
c
)
>
0 and ℝ0>1
(b)
H
(
I
*
,
c
)
<
0 and ℝ0>1
(c)
H
(
I
*
,
c
)
>
0 and ℝ0<1
(d)
H
(
I
*
,
c
)
<
0 and ℝ0<1From Theorem5, we know that the positive equilibrium E* of system (6) is a center-type nonhyperbolic equilibrium when H(I*,c)=0. Hence, system (6) may undergo Hopf bifurcation in this case. To determine the stability of the endemic equilibrium and direction of Hopf bifurcation in this case, we must compute the Lyapunov coefficients of the equilibrium. We first translate the endemic equilibrium E* of system (6) to the origin with x=S-S*,y=I-I*. Then system (6) in a neighborhood U of the origin can be written as
(28)dxdt=a10x+a01y+a11xy,dydt=b10x+b01y+b11xy+b02y2+b03y3+O(y4),
where a10=-d-βI*,a01=-βS*,a11=-β,b10=βI*,b01=cI*/(b+I*)2,b11=β,b02=bc/(b+I*)3,b03=-bc/(b+I*)4 and O(y4) is the same order infinity. Hence, using the formula of the Lyapunov number σ for the focus at the origin of (28) in [13], we have
(29)σ=-3π2a01Δ3/2×[a10b10(a112+a11b02)+a10a01(b112+a11b02)-2a10b10b022+(a01b10-2a102)b11b02-3b10b03(a102+a01b10)],
where Δ=a10b01-a01b10=I*G(I*,c)>0.By numerical simulation, if we fixd=0.02,ν=26, we know there exist parameter values (A,b,c,β)=(1050,500,128228,0.5), which satisfies ℝ0=92.93>1 and H(I*,c)=0 such that σ=-97.09<0. On the other hand, there exist values (A,b,c,β)=(1050,0.1,357.3097,0.5), which satisfy ℝ0=7.29>1 and H(I*,c)=0 too, but σ=0.0015>0. Therefore, there exists an open set V1 in the parameter space (A,b,c,β), such that H(I*,c)=0,ℝ0>1 and σ<0; that is,
(30)V1={(A,b,c,β):A>0,b>0,c>0,β>0,ℝ0>1,H(I*,c)=0,σ<0}.
And there exists another open set V2 in the parameter space (A,b,c,β), such that H(I*,c)=0,ℝ0>1 and σ>0; that is,
(31)V2={(A,b,c,β):A>0,b>0,c>0,β>0,ℝ0>1,H(I*,c)=0,σ>0}.
By [13], we have the following theorem.Theorem 6.
Supposeℝ0>1 and there exists ck>0 such that H(I*,ck)=0 and (∂H(I*,c)/∂c)|c=ck≠0 for the endemic equilibrium E* of the model. If σ≠0, then a curve of periodic solutions bifurcates from E* such that the following happens. (1) Suppose (A,b,ck,β)∈V1; the model undergoes a supercritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck<0 and backward supercritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck>0. (2) Suppose (A,b,ck,β)∈V2; the model undergoes a subcritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck<0 and backward subcritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck>0.


### 4.2. Case 2:ℝ^0<ℝ0<1

Whenℝ^0<ℝ0<1, model (6) has two endemic equilibria E*(S*,I*) and E*(S*,I*) when a1<0(A>A1(c)). Their coordinates satisfy (7) and I*<I*. In this case, we will discuss the stability of endemic equilibrium, Hopf bifurcation, similarly as we did in the previous subsection. In addition, backward bifurcation will be discussed too.It is easy to see that the Jacobian matrix (6) at E* and E* are the same as (24) if one replaces I* with I*. We have the following result.Theorem 7.
Supposeℝ^0<ℝ0<1 and a1<0. Then the endemic equilibrium E* of (6) is a saddle; E* is a stable node or focus when H(I*,c)>0;E* is an unstable node or focus when H(I*,c)<0;E* is a center of the linear system when H(I*,c)=0.Proof.
Becausedg(I*)/dI-dh(I*)/dI<0,dg(I*)/dI-dh(I*)/dI>0, we have G(I*,c)<0,G(I*,c)>0 from (27). Then E* is a saddle, E* is a stable node or focus when H(I*,c)>0,E* is an unstable node or focus when H(I*,c)<0,andE* is a center of the linear system when H(I*,c)=0.Here we give some numerical simulations to show thatH(I*,c) may be positive or negative for different parameters. For example, we fix A=15000,β=0.0005,d=0.02,b=0.01,ν=34. When c=7.16, we have ℝ0=0.50. System (6) has two endemic equilibria, E*(S*,I*)=(749793.72,0.011)andE*(S*,I*)=(68075.74,400.69). We get H(I*,c)=0.20247. From Figure 3(c), E* is locally stable. When c=201, we have ℝ0=0.019. System (6) has two endemic equilibria too, E*(S*,I*)=(739115.82,0.589),andE*(S*,I*)=(69059.22,394.41). In this case, H(I*,c)=-0.29239. From Figure 3(d), we can see that E* is unstable.From Theorem7, we know that the endemic equilibrium E* of system (6) is a center-type nonhyperbolic equilibrium when H(I*,c)=0. Hence, system (6) may undergo Hopf bifurcation in this case. To determine the stability of the endemic equilibrium and direction of Hopf bifurcation in this case, we must compute the Lyapunov coefficients σ of the equilibrium E* as we did in Case 1.Denote nonempty setV3 and V4 as
(32)V3={ℝ^0(A,b,c,β):A>0,b>0,c>0,β>0,ℝ^0(c)<ℝ0<1,H(I*,c)=0,σ<0},(33)V4={ℝ^0(A,b,c,β):A>0,b>0,c>0,β>0,ℝ^0(c)<ℝ0<1,H(I*,c)=0,σ>0}.Theorem 8.
Supposeℝ^0<ℝ0<1,a1<0(A>A1(c)) and there exists ck>0 such that H(I*,ck)=0 and (∂H(I*,c)/∂c)|c=ck≠0 for the endemic equilibrium E* of model (6). If σ≠0, then a curve of periodic solutions bifurcates from E2 such that the following happens. (1) Suppose (A,b,ck,β)∈V3; then the model undergoes a supercritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck<0 and backward supercritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck>0. (2) Suppose (A,b,ck,β)∈V4; then the model undergoes a subcritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck<0 and backward subcritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck>0.According to the aforementioned results, we know that whend-bβ>0 and c>c*, besides the basic reproduction number ℝ0, there exists a subthreshold condition for model (3). If ℝ0≥1, model (3) has a unique endemic equilibrium which is a node or focus. If ℝ^0<ℝ0<1, model (3) has two endemic equilibria; one is a node or focus and the other is a saddle point. If ℝ0=ℝ^0, model (3) has a unique endemic equilibria of multiplicity two. Otherwise, there is no endemic equilibrium. This situation corresponds to a backward bifurcation (Figure 4). We translate the results into the following theorem.Figure 4
A backward bifurcation atℝ0=1 with a dashed curve for the location of the saddle point and solid curve for the location of the other endemic equilibria.Theorem 9.
Ifd-bβ>0 and c>c*, model (3) undergoes a backward bifurcation at ℝ0=1.Remark 10.
Since the existence of two equilibria and limit cycle under some conditions, model (3) may exhibit two important mathematical phenomena, bistability and periodicity, which can be observed in data for some infectious diseases in history.


## 4.1. Case 1:ℝ0>1

Whenℝ0>1, the model (6) has a unique endemic equilibrium E*. From (24) to (26), we have the following.Theorem 5.
Supposeℝ0>1. Then the disease endemic equilibrium E* of (6) is a stable node or focus when H(I*,c)>0;E* is an unstable node or focus when H(I*,c)<0, and (6) has at least one closed orbit in Ω;E* is a center of the linear system when H(I*,c)=0.Proof.
RewriteG(I,c) as
(27)G(I,c)=β(d+βI)[-cβ(b+I)2+βA(d+βI)2]=β(d+βI)[dg(I)dI-dh(I)dI].
Note that dg(I*)/dI-dh(I*)/dI>0; hence G(I*,c)>0. The unique endemic equilibrium E* is a stable node or focus when H(I*,c)>0,E* is an unstable node or focus when H(I*,c)<0;E* is a center of the linear system when H(I*,c)=0. If H(I*,c)<0, then E is an unstable node or focus. System (6) has at least one closed orbit in Ω from the Poincaré-Bendixson Theorem.As an example, we fixA=15000,β=0.0005,d=0.02,b=1,ν=34. When c=30, then ℝ0=5.86; the model (6) has unique endemic equilibrium E(S*,I*)=(68189.645,399.949) and H(I*,c)=0.14534>0. From Figure 3(a), we can see that the endemic equilibrium E* is stable. When c=126.99, then ℝ0=2.33; the model also has a unique endemic equilibrium E(S*,I*)=(68678.432,396.818) but H(I*,c)=-0.01<0. From Figure 3(b), we can see that the endemic equilibrium E* is unstable and there must be a stable limit cycle around E(S*,I*).The phase portraits of the model where the directions of trajectories are counterclockwise.
(a)
H
(
I
*
,
c
)
>
0 and ℝ0>1
(b)
H
(
I
*
,
c
)
<
0 and ℝ0>1
(c)
H
(
I
*
,
c
)
>
0 and ℝ0<1
(d)
H
(
I
*
,
c
)
<
0 and ℝ0<1From Theorem5, we know that the positive equilibrium E* of system (6) is a center-type nonhyperbolic equilibrium when H(I*,c)=0. Hence, system (6) may undergo Hopf bifurcation in this case. To determine the stability of the endemic equilibrium and direction of Hopf bifurcation in this case, we must compute the Lyapunov coefficients of the equilibrium. We first translate the endemic equilibrium E* of system (6) to the origin with x=S-S*,y=I-I*. Then system (6) in a neighborhood U of the origin can be written as
(28)dxdt=a10x+a01y+a11xy,dydt=b10x+b01y+b11xy+b02y2+b03y3+O(y4),
where a10=-d-βI*,a01=-βS*,a11=-β,b10=βI*,b01=cI*/(b+I*)2,b11=β,b02=bc/(b+I*)3,b03=-bc/(b+I*)4 and O(y4) is the same order infinity. Hence, using the formula of the Lyapunov number σ for the focus at the origin of (28) in [13], we have
(29)σ=-3π2a01Δ3/2×[a10b10(a112+a11b02)+a10a01(b112+a11b02)-2a10b10b022+(a01b10-2a102)b11b02-3b10b03(a102+a01b10)],
where Δ=a10b01-a01b10=I*G(I*,c)>0.By numerical simulation, if we fixd=0.02,ν=26, we know there exist parameter values (A,b,c,β)=(1050,500,128228,0.5), which satisfies ℝ0=92.93>1 and H(I*,c)=0 such that σ=-97.09<0. On the other hand, there exist values (A,b,c,β)=(1050,0.1,357.3097,0.5), which satisfy ℝ0=7.29>1 and H(I*,c)=0 too, but σ=0.0015>0. Therefore, there exists an open set V1 in the parameter space (A,b,c,β), such that H(I*,c)=0,ℝ0>1 and σ<0; that is,
(30)V1={(A,b,c,β):A>0,b>0,c>0,β>0,ℝ0>1,H(I*,c)=0,σ<0}.
And there exists another open set V2 in the parameter space (A,b,c,β), such that H(I*,c)=0,ℝ0>1 and σ>0; that is,
(31)V2={(A,b,c,β):A>0,b>0,c>0,β>0,ℝ0>1,H(I*,c)=0,σ>0}.
By [13], we have the following theorem.Theorem 6.
Supposeℝ0>1 and there exists ck>0 such that H(I*,ck)=0 and (∂H(I*,c)/∂c)|c=ck≠0 for the endemic equilibrium E* of the model. If σ≠0, then a curve of periodic solutions bifurcates from E* such that the following happens. (1) Suppose (A,b,ck,β)∈V1; the model undergoes a supercritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck<0 and backward supercritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck>0. (2) Suppose (A,b,ck,β)∈V2; the model undergoes a subcritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck<0 and backward subcritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck>0.


## 4.2. Case 2:ℝ^0<ℝ0<1

Whenℝ^0<ℝ0<1, model (6) has two endemic equilibria E*(S*,I*) and E*(S*,I*) when a1<0(A>A1(c)). Their coordinates satisfy (7) and I*<I*. In this case, we will discuss the stability of endemic equilibrium, Hopf bifurcation, similarly as we did in the previous subsection. In addition, backward bifurcation will be discussed too.It is easy to see that the Jacobian matrix (6) at E* and E* are the same as (24) if one replaces I* with I*. We have the following result.Theorem 7.
Supposeℝ^0<ℝ0<1 and a1<0. Then the endemic equilibrium E* of (6) is a saddle; E* is a stable node or focus when H(I*,c)>0;E* is an unstable node or focus when H(I*,c)<0;E* is a center of the linear system when H(I*,c)=0.Proof.
Becausedg(I*)/dI-dh(I*)/dI<0,dg(I*)/dI-dh(I*)/dI>0, we have G(I*,c)<0,G(I*,c)>0 from (27). Then E* is a saddle, E* is a stable node or focus when H(I*,c)>0,E* is an unstable node or focus when H(I*,c)<0,andE* is a center of the linear system when H(I*,c)=0.Here we give some numerical simulations to show thatH(I*,c) may be positive or negative for different parameters. For example, we fix A=15000,β=0.0005,d=0.02,b=0.01,ν=34. When c=7.16, we have ℝ0=0.50. System (6) has two endemic equilibria, E*(S*,I*)=(749793.72,0.011)andE*(S*,I*)=(68075.74,400.69). We get H(I*,c)=0.20247. From Figure 3(c), E* is locally stable. When c=201, we have ℝ0=0.019. System (6) has two endemic equilibria too, E*(S*,I*)=(739115.82,0.589),andE*(S*,I*)=(69059.22,394.41). In this case, H(I*,c)=-0.29239. From Figure 3(d), we can see that E* is unstable.From Theorem7, we know that the endemic equilibrium E* of system (6) is a center-type nonhyperbolic equilibrium when H(I*,c)=0. Hence, system (6) may undergo Hopf bifurcation in this case. To determine the stability of the endemic equilibrium and direction of Hopf bifurcation in this case, we must compute the Lyapunov coefficients σ of the equilibrium E* as we did in Case 1.Denote nonempty setV3 and V4 as
(32)V3={ℝ^0(A,b,c,β):A>0,b>0,c>0,β>0,ℝ^0(c)<ℝ0<1,H(I*,c)=0,σ<0},(33)V4={ℝ^0(A,b,c,β):A>0,b>0,c>0,β>0,ℝ^0(c)<ℝ0<1,H(I*,c)=0,σ>0}.Theorem 8.
Supposeℝ^0<ℝ0<1,a1<0(A>A1(c)) and there exists ck>0 such that H(I*,ck)=0 and (∂H(I*,c)/∂c)|c=ck≠0 for the endemic equilibrium E* of model (6). If σ≠0, then a curve of periodic solutions bifurcates from E2 such that the following happens. (1) Suppose (A,b,ck,β)∈V3; then the model undergoes a supercritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck<0 and backward supercritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck>0. (2) Suppose (A,b,ck,β)∈V4; then the model undergoes a subcritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck<0 and backward subcritical Hopf bifurcation if (∂H(I*,c)/∂c)|c=ck>0.According to the aforementioned results, we know that whend-bβ>0 and c>c*, besides the basic reproduction number ℝ0, there exists a subthreshold condition for model (3). If ℝ0≥1, model (3) has a unique endemic equilibrium which is a node or focus. If ℝ^0<ℝ0<1, model (3) has two endemic equilibria; one is a node or focus and the other is a saddle point. If ℝ0=ℝ^0, model (3) has a unique endemic equilibria of multiplicity two. Otherwise, there is no endemic equilibrium. This situation corresponds to a backward bifurcation (Figure 4). We translate the results into the following theorem.Figure 4
A backward bifurcation atℝ0=1 with a dashed curve for the location of the saddle point and solid curve for the location of the other endemic equilibria.Theorem 9.
Ifd-bβ>0 and c>c*, model (3) undergoes a backward bifurcation at ℝ0=1.Remark 10.
Since the existence of two equilibria and limit cycle under some conditions, model (3) may exhibit two important mathematical phenomena, bistability and periodicity, which can be observed in data for some infectious diseases in history.


## 5. Discussion

In this paper, we proposed a SIR epidemic model with saturation recovery to understand the impact of limited medical resource on infectious disease transmission. Existence of equilibria was obtained under different conditions and their stabilities were analyzed. Bifurcations, including backward bifurcation and Hopf bifurcation, were analyzed too.The basic reproduction numberℝ0, which gives the expected number of new infections from one infectious individual over the duration of the infectious period, given that all other members of the population are susceptible, in terms of the model parameters, was identified. Our results suggested that for ℝ0>1, the disease-free equilibrium point was unstable while the endemic equilibrium emerged as a unique equilibrium point which implies that disease never dies out epidemiologically. Therefore, bringing the basic reproduction number below 1 is essential. Nevertheless, if ℝ0<1, from Theorems 2, 3, and 9, the disease-free equilibrium point was not always globally stable and the model may undergo a backward bifurcation, bistability, and periodicity, which implies that reducing the basic reproduction number below 1 is not enough to control the disease in some cases. In order to eradicate the disease, we have a subthreshold number ℝ^0 and it is now necessary to reduce ℝ0 to a value less than ℝ^0.c, the maximum treatment per unit of time, is related to the maximum treatment capacity in a city or region. Note that ℝ0 is a monotone decreasing function of c. Therefore, under the condition of d-bβ≤0, increasing c to a value such that ℝ0<1 is sufficient to eliminate the disease. While under the condition of d-bβ>0, increasing c to a value such that ℝ0<ℝ^0 is necessary. Decreasing c to a value at which ℝ^0<ℝ0<1, then whether the disease will break out may depend on the initial conditions. There is a region such that the disease dies out if the initial position lies in this region; otherwise, it tends to an endemic equilibrium E* or a periodic orbit around E*. Since the eventual behavior is related to the initial position, this model may be realistic and useful and we should pay more attention to the initial states of the disease (see Figure 5). If c is small enough such that ℝ0>1, the disease must break out. Therefore, a sufficient number of sickbeds and other medical resources are very important for the disease control and eradication.Different initial states bring different results. The parameters areA=15000,β=0.0005,d=0.02,b=0.01,ν=34,c=7.16, and ℝ0=0.50. The unit of time t is day.
(a)
The time course ofI(t) when S(0)=68075,I(0)=400
(b)
The time course ofI(t) when S(0)=68075,I(0)=3The mathematical analysis of the SIR model has highlighted two important mathematical phenomena observed in data for some infectious diseases: bistability and periodicity. In the bistability phenomenon, the model exhibits multiple endemic equilibria even whenℝ0<1. A clinical study of measles during an endemic in Poland shows that despite high vaccination coverage needed for eradiation since the 1980s, an epidemic of measles with 2255 reported cases occurred between November 1997 and 1998 ([1, 14]). The results of this study confirm that reducing ℝ0 to values less than unity may fail to control the disease. Secondly, the model undergoes oscillatory behavior under certain conditions. The existence of limit cycles confirms such behavior which has been reported in many studies on the dynamics of some infectious diseases such as measles, rubella, and so forth. ([1, 10, 15]).There are many interesting problems related to the models with situation recovery. Some of them for the model will be studied in future. Important dynamical systems features of our models, for example, saddle-node bifurcation, and degenerate Hopf bifurcations, should be discussed.



---

*Source: 314958-2013-04-30.xml*