Rich Dynamics of an Epidemic Model with Saturation Recovery

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 ℝ₀ 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.

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 Section 2, 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) to model the treated part which is increasing for small infectives and approaches the maximum for large infectives. Hence 𝑐, as the limit of ℎ(𝐼) as 𝐼 tends to infinity, is the maximum of the treatment capacity in a region, and 𝑏, the infected size at which is 50% saturation (ℎ(𝑏)=𝑐/2), measures how soon saturation occurs. Here ℎ(𝐼) satisfies ℎ(0)=0,ℎ(∞)=𝑐 and 𝑑ℎ/𝑑𝐼>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. Let 𝑆(𝑡),𝐼(𝑡), and 𝑅(𝑡) denote the number of susceptible, infective, and recovered individuals at time 𝑡, 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) where all the parameters are positive and(i)𝐴is the recruitment rate of susceptible population;(ii)𝑑 is natural death rate and 𝜈 is the disease-induced death rate;(iii)𝑐 is the maximum of treatment per unit of time, and 𝑏, the infected size at which is 50% saturation (ℎ(𝑏)=𝑐/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)

Model (3) has one disease-free equilibrium at 𝐸₀=(𝐴/𝑑,0). Using the formulae in [12], a straightforward calculation gives the reproduction number: ℝ₀=𝐴𝑏𝛽𝑑().𝑏𝑑+𝑏𝜈+𝑐 (4)

The disease-free equilibrium 𝐸₀ has two eigenvalues −𝑑 and ℝ₀−1. Therefore we have the following proposition.

Proposition 1. For the model (3), the disease-free equilibrium 𝐸₀ is locally asymptotically stable if ℝ₀<1 and unstable if ℝ₀>1.

3. Existence of the Endemic Equilibrium

Let ℎ(𝐼)=𝐴(𝐼)=𝑑+𝛽𝐼,𝑔𝑑+𝜈𝛽+𝑐𝛽().𝑏+𝐼 (5) Then we can rewrite the model (3) as 𝑑𝑆((𝐼,𝑑𝑡=−𝑑+𝛽𝐼𝑆−ℎ))𝑑𝐼(𝐼.𝑑𝑡=𝛽𝐼𝑆−𝑔) (6)

Let the right hand side of (6) be zero. If an endemic equilibrium exists, its (𝑆,𝐼) coordinates must satisfy 𝑆=𝑔(𝐼),𝑆=ℎ(𝐼). (7) We note that lim_𝐼→∞ℎ(𝐼)=0,lim_𝐼→∞𝑔(𝐼)=(𝑑+𝜈)/𝛽,𝑑ℎ/𝑑𝐼<0,𝑑𝑔/𝑑𝐼<0, and 𝑑²ℎ/𝑑𝐼²>0,𝑑²𝑔/𝑑𝐼²>0. In addition, ℝ₀=ℎ(0)/𝑔(0).

Eliminating 𝑆 from (7) gives an equation of the form 𝐼²+𝑎₁𝐼+𝑎₂=0, (8) where 𝑎₁=𝑏𝛽+𝑑(𝑑+𝜈)+𝑐𝛽−𝐴𝛽𝛽(),𝑎𝑑+𝜈₂=𝑏𝑑(𝑑+𝜈)+𝑐𝑑−𝛽𝐴𝑏𝛽()=()𝑑+𝜈𝑏𝑑𝑑+𝜈+𝑐𝑑𝛽(𝑑+𝜈)1−ℝ₀. (9)

If ℝ₀>1, then 𝑎₂<0, and there is a unique positive root for (8) which implies that a unique endemic equilibrium 𝐸^∗(𝑆^∗,𝐼^∗) exists (see Figure 1(b)).

(a) No positive equilibrium

(b) Unique positive equilibrium

(c) Unique equilibrium of multiplicity 2

(d) Two positive equilibria

Figure 1 

Existence and number of endemic equilibria.

If ℝ₀=1, then 𝑎₂=0 and there is a unique nonzero solution of (8) 𝐼=−𝑎₁ which is positive if and only if 𝑎₁<0. Then, there is a unique endemic equilibrium 𝐸^∗(𝑆^∗,𝐼^∗) when 𝑎₁<0, and there are not endemic equilibria when 𝑎₁≥0.

If ℝ₀<1, then 𝑎₂>0. For (8) to have at least one positive root we must have 𝑎₁<0,Δ≜𝑎²₁−4𝑎₂⩾0. (10) Solving Δ=0 in terms of ℝ₀, one gets ℝ₀=ℝ₀ where ℝ₀(𝑐)=𝐴𝑏(𝑑+𝜈)𝑏(𝑑+𝜈)+√𝐴+√𝑐²×((𝑏(𝑏()(+()𝑐+𝑏𝑑+𝜈))𝑑+𝜈)(𝑑+𝜈+2𝐴+𝑐))𝐴−𝑐²⁻¹. (11) One can verify that, provided 𝑎₁<0, model (3) has exactly 0, 1, and 2 endemic equilibria as shown in Figures 1(a), 1(c), and 1(d) for ℝ₀<ℝ₀,ℝ₀=ℝ₀,andℝ₀>ℝ₀, respectively.

In summary, regarding the existence of endemic equilibrium, we have the following.

Theorem 2. For model (3), one has the following.(1)If ℝ₀>1, there exists a unique positive equilibrium 𝐸^∗(𝑆^∗,𝐼^∗).(2)If ℝ₀=1; there is a positive equilibrium 𝐸^∗(𝑆^∗,𝐼^∗) when 𝑎₁<0; otherwise there is no positive equilibrium.(3)If ℝ₀<1 and 𝑎₁≥0, there is no positive equilibrium.(4)If ℝ₀<ℝ₀<1 and 𝑎₁<0, there are two positive equilibria 𝐸^∗ and 𝐸_∗.(5)If ℝ₀=ℝ₀ and 𝑎₁<0,𝐸^∗ and 𝐸_∗ coalesce together as a unique equilibrium of multiplicity two.(6)If ℝ₀<ℝ₀ and 𝑎₁<0, there is no positive equilibrium.
When exist, 𝐸^∗(𝑆^∗,𝐼^∗) and 𝐸_∗(𝑆_∗,𝐼_∗) are the corresponding equilibria, and 𝑆^∗=ℎ𝐼^∗,𝑆_∗=ℎ𝐼_∗,𝐼^∗=−𝑎₁+√Δ2,𝐼_∗=−𝑎₁−√Δ2. (12)

𝑐, the maximum treatment per unit of time, is related to the treatment capacity in a region and the recruitment rate of susceptible individuals 𝐴 can be changed by vaccination, and so forth; they are all important for disease control. In the following section, we choose 𝑐 and 𝐴 as parameters to discuss the distribution of the equilibria of system (6) in plane (𝑐,𝐴). Solving equations 𝑎₁=0 and 𝑎₂=0 for 𝐴, we can get 𝐴=𝐴₁(𝑐) and 𝐴=𝐴₂(𝑐), respectively, where 𝐴₁(𝑐)𝑑=𝑐+𝑏+𝛽(),𝐴𝑑+𝜈₂(𝑐)=𝑑𝑑()𝑏𝛽𝑐+𝑑+𝜈𝛽. (13) Then 𝑎₁=𝐴₁(𝑐)1−𝐴,𝑎𝑑+𝜈₂=𝐴₂(𝑐)𝑏−𝐴.𝑑+𝜈 (14)

Obviously, 𝑎₁>0 (or 𝑎₁<0) when 𝐴<𝐴₁(𝑐) (or 𝐴>𝐴₁(𝑐)) and 𝑎₂>0 (or 𝑎₂<0) when 𝐴<𝐴₂(𝑐) (or 𝐴>𝐴₂(𝑐)).

It is easy to see that if 𝐴>𝐴₂(𝑐), (7) has a unique positive solution.

Now we consider the case when 𝐴⩽𝐴₂(𝑐) (the case 𝑎₂≥0). Writing the discriminant Δ=𝑎²₁−4𝑎₂ as a function of 𝑐 and 𝐴, we have Δ()=)𝑐,𝐴𝑏𝛽+𝑑(𝑑+𝜈+𝑐𝛽−𝐴𝛽𝛽(𝑑+𝜈)²()−4𝑏𝑑𝑑+𝜈+𝑐𝑑−𝛽𝐴𝑏𝛽().𝑑+𝜈 (15) If Δ=0 defines curves in the first quadrant of the (𝑐,𝐴) plane, then it must satisfy 𝐴⩽𝐴₂(𝑐). Furthermore, if (8) has nonnegative roots, we need 𝐴₁(𝑐)<𝐴⩽𝐴₂(𝑐) and Δ⩾0. Calculate the difference of 𝐴₂(𝑐) and 𝐴₁(𝑐): 𝐴₂(𝑐)−𝐴₁(𝑐)=𝑑−𝑏𝛽().𝑏𝛽𝑐−𝑏𝑑+𝜈 (16)

If 𝑑−𝑏𝛽⩽0, then 𝐴₂(𝑐)<𝐴₁(𝑐) and system (6) does not have any endemic equilibrium when 𝐴⩽𝐴₂(𝑐).

If 𝑑−𝑏𝛽>0, one easily gets that 𝑐=𝑐^∗ from 𝐴₂(𝑐)−𝐴₁(𝑐)=0, where 𝑐^∗=𝑏²𝛽(𝑑+𝜈).𝑑−𝑏𝛽 (17) We have 𝑐⩽𝑐^∗⟺𝐴₁(𝑐)⩾𝐴₂(𝑐),𝑐>𝑐^∗⟺𝐴₁(𝑐)<𝐴₂(𝑐). (18) Therefore, when 𝑐⩽𝑐^∗, system (6) does not have any endemic equilibrium in this case. If 𝑐>𝑐^∗, then system (6) has two equilibria when 𝐴₁(𝑐)<𝐴<𝐴₂(𝑐) and Δ>0; system (6) has a unique endemic equilibrium when 𝐴=𝐴₂(𝑐); system (6) has a unique endemic equilibrium of multiplicity 2 when 𝐴₁(𝑐)<𝐴<𝐴₂(𝑐) and Δ=0; and system (6) has no equilibrium when Δ<0.

Solving the equation Δ=0 for 𝐴, one gets two curves Γ₁∶𝐴=𝐴₃(𝑐),Γ₂∶𝐴=𝐴₄(𝑐), (19) in the (𝑐,𝐴) plane, where 𝐴₃(𝑐)=𝑑−𝑏𝛽(𝑑+𝜈)+𝑐𝛽+2𝑐𝛽𝑑−𝑏𝛽(𝑑+𝜈)𝛽,𝐴₄(𝑐)=𝑑−𝑏𝛽(𝑑+𝑣)+𝑐𝛽−2𝑐𝛽𝑑−𝑏𝛽(𝑑+𝜈)𝛽. (20) For any point of (𝑐,𝐴)∈Γ₂, we have 𝐴<𝐴₁(𝑐) by simple calculation. The two lines 𝐴=𝐴₁(𝑐),𝐴=𝐴₂(𝑐) joint the curve Γ₁∶𝐴=𝐴₃(𝑐) at a point 𝑃(𝐴₁(𝑐^∗),𝑐^∗) when 𝑐=𝑐^∗. For any 𝑐>𝑐^∗, the curve Γ₁ lies above 𝐴=𝐴₁(𝑐) and under 𝐴=𝐴₂(𝑐) as shown in Figure 2. In fact, we have Δ<0 for the point (𝑐,𝐴) on 𝐴=𝐴₁(𝑐) and Δ>0 for the point (𝑐,𝐴) on 𝐴=𝐴₂(𝑐), and hence the curve Γ₁ lies between the two lines 𝐴=𝐴₁(𝑐) and 𝐴=𝐴₂(𝑐) when 𝑐>𝑐^∗. As is shown in Figure 2, in the case when 𝑑−𝑏𝛽>0, the curve 𝐴=𝐴₂(𝑐) and the curve segment Γ₁ subdivide the positive cone of the parameters plane (𝑐,𝐴) into three regions 𝐷₀,𝐷₁ and 𝐷₂, where 𝐷₀=(𝑐,𝐴)∣0<𝑐⩽𝑐^∗,0<𝐴⩽𝐴₂(𝑐∪))(𝑐,𝐴∣𝑐>𝑐^∗,0<𝐴<𝐴₃(𝑐,𝐷₁)=(𝑐,𝐴)∣0<𝑐⩽𝑐^∗,𝐴>𝐴₂(𝑐∪))(𝑐,𝐴∣𝑐>𝑐^∗,𝐴⩾𝐴₂(𝑐,𝐷₂)=(𝑐,𝐴)∣𝑐>𝑐^∗,𝐴₃(𝑐)<𝐴<𝐴₂(𝑐). (21) For parameter (𝑐,𝐴) 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 (𝑐,𝐴) when 𝑑−𝑏𝛽>0. There are 0, 1, 2 equilibria in 𝐷0,𝐷1,𝐷2, respectively.

Theorem 3. For the model (6), with 𝑐^∗,𝐴₁(𝑐),𝐴₂(𝑐), and 𝐴₃(𝑐) defined as previously mentioned,(1)if 𝑑−𝑏𝛽>0, as shown in Figure 2, one has the following.(a) For (𝑐,𝐴)∈𝐷₀, the model (6) has no equilibrium.(b) For (𝑐,𝐴)∈𝐷₁, the model (6) has a unique equilibrium 𝐸^∗.(c) For (𝑐,𝐴)∈𝐷₂, the model (6) has two equilibria 𝐸_∗ and 𝐸^∗.(d) For (𝑐,𝐴)∈Γ₁, when 𝑐>𝑐^∗,𝐸_∗ and 𝐸^∗ coalesce at a unique endemic equilibrium of multiplicity 2.(e) At the point 𝑃(𝐴₁(𝑐^∗),𝑐^∗), the model (6) has no endemic equilibrium; the disease-free equilibrium is of multiplicity 2.(2)On the other hand, if 𝑑−𝑏𝛽⩽0, then(a) if 𝐴>𝐴₂(𝑐), there is a unique endemic equilibrium 𝐸^∗;(b) if 𝐴⩽𝐴₂(𝑐), there is no endemic equilibrium.
When exist, 𝐸^∗(ℎ(𝐼^∗),𝐼^∗) and 𝐸_∗(ℎ(𝐼_∗),𝐼_∗) are the corresponding equilibria, and 𝐼^∗=12−𝑎₁+Δ(,𝐼𝑐,𝐴)_∗=12−𝑎₁−Δ(𝑐,𝐴). (22)

Now we can prove the global stability of the disease-free equilibrium 𝐸₀ when (𝑐,𝐴)∈𝐷₀ according to Theorem 3.

Corollary 4. If (𝑐,𝐴)∈𝐷₀, then the disease-free equilibrium 𝐸₀ of (6) is globally asymptotically stable.

Proof. Let Ω={(𝑆,𝐼)∣𝑆,𝐼≥0,𝑆+𝐼≤𝐴/𝑑}. We will prove that Ω is a positively invariant set of (6). By (6), we have (𝑑𝑆/𝑑𝑡)|_𝑆=0=𝐴>0 and (𝑑𝐼/𝑑𝑡)|_𝐼=0=0. Denote 𝑁(𝑡)=𝑆(𝑡)+𝐼(𝑡); then 𝑑𝑁𝑑𝑡_{𝑁=𝐴/𝑑}≤𝐴−𝑑𝑁|_{𝑁=𝐴/𝑑}=0. (23) Hence Ω is a positively invariant set of (6), and it attracts all the positive orbits of (6) state in 𝑅⁺₂. We have shown that (6) does not have any positive equilibria when (𝑐,𝐴)∈𝐷₀. 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 𝐸₀ is the unique equilibrium in Ω, the local stability of 𝐸₀ implies that the 𝜔-limit set of every solution with an initial point in 𝑅⁺₂ is 𝐸₀. Hence, the disease-free equilibrium 𝐸₀ 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 𝐸^∗(𝑆^∗,𝐼^∗) gives 𝐽=−𝑑−𝛽𝐼^∗−𝛽𝑆^∗𝛽𝐼^∗𝑐𝐼^∗(𝑏+𝐼^∗)². (24) Then the characteristic equation about 𝐸^∗ is given by 𝜆²+𝐻𝐼^∗,𝑐𝜆+𝐼^∗𝐺𝐼^∗,𝑐=0, (25) where 𝐻𝐼^∗,𝑐=𝑑+𝛽𝐼^∗−𝑐𝐼^∗(𝑏+𝐼^∗)²,𝐺(𝐼,𝑐)=𝐴𝛽²𝑑+𝛽𝐼^∗−𝑐𝑑+𝛽𝐼^∗(𝑏+𝐼^∗)². (26)

4.1. Case 1: ℝ₀>1

When ℝ₀>1, the model (6) has a unique endemic equilibrium 𝐸^∗. From (24) to (26), we have the following.

Theorem 5. Suppose ℝ₀>1. Then the disease endemic equilibrium 𝐸^∗ of (6) is a stable node or focus when 𝐻(𝐼^∗,𝑐)>0;𝐸^∗ is an unstable node or focus when 𝐻(𝐼^∗,𝑐)<0, and (6) has at least one closed orbit in Ω;𝐸^∗ is a center of the linear system when 𝐻(𝐼^∗,𝑐)=0.

Proof. Rewrite 𝐺(𝐼,𝑐) as 𝐺()−𝑐𝐼,𝑐=𝛽𝑑+𝛽𝐼𝛽(𝑏+𝐼)²+𝛽𝐴𝑑+𝛽𝐼²(𝐼)=𝛽𝑑+𝛽𝐼𝑑𝑔−(𝐼)𝑑𝐼𝑑ℎ.𝑑𝐼 (27) Note that 𝑑𝑔(𝐼^∗)/𝑑𝐼−𝑑ℎ(𝐼^∗)/𝑑𝐼>0; hence 𝐺(𝐼^∗,𝑐)>0. The unique endemic equilibrium 𝐸^∗ is a stable node or focus when 𝐻(𝐼^∗,𝑐)>0,𝐸^∗ is an unstable node or focus when 𝐻(𝐼^∗,𝑐)<0;𝐸^∗ is a center of the linear system when 𝐻(𝐼^∗,𝑐)=0. If 𝐻(𝐼^∗,𝑐)<0, then 𝐸 is an unstable node or focus. System (6) has at least one closed orbit in Ω from the Poincaré-Bendixson Theorem.

As an example, we fix 𝐴=15000,𝛽=0.0005,𝑑=0.02,𝑏=1,𝜈=34. When 𝑐=30, then ℝ₀=5.86; the model (6) has unique endemic equilibrium 𝐸(𝑆^∗,𝐼^∗)=(68189.645,399.949) and 𝐻(𝐼^∗,𝑐)=0.14534>0. From Figure 3(a), we can see that the endemic equilibrium 𝐸^∗ is stable. When 𝑐=126.99, then ℝ₀=2.33; the model also has a unique endemic equilibrium 𝐸(𝑆^∗,𝐼^∗)=(68678.432,396.818) but 𝐻(𝐼^∗,𝑐)=−0.01<0. From Figure 3(b), we can see that the endemic equilibrium 𝐸^∗ is unstable and there must be a stable limit cycle around 𝐸(𝑆^∗,𝐼^∗).

(a) 𝐻(𝐼∗,𝑐)>0 and ℝ0>1

(b) 𝐻(𝐼∗,𝑐)<0 and ℝ0>1

(c) 𝐻(𝐼∗,𝑐)>0 and ℝ0<1

(d) 𝐻(𝐼∗,𝑐)<0 and ℝ0<1

Figure 3 

The phase portraits of the model where the directions of trajectories are counterclockwise.

From Theorem 5, we know that the positive equilibrium 𝐸^∗ of system (6) is a center-type nonhyperbolic equilibrium when 𝐻(𝐼^∗,𝑐)=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 𝐸^∗ of system (6) to the origin with 𝑥=𝑆−𝑆^∗,𝑦=𝐼−𝐼^∗. Then system (6) in a neighborhood 𝑈 of the origin can be written as 𝑑𝑥𝑑𝑡=𝑎₁₀𝑥+𝑎₀₁𝑦+𝑎₁₁𝑥𝑦,𝑑𝑦𝑑𝑡=𝑏₁₀𝑥+𝑏₀₁𝑦+𝑏₁₁𝑥𝑦+𝑏₀₂𝑦²+𝑏₀₃𝑦³+𝑂𝑦⁴, (28) where 𝑎₁₀=−𝑑−𝛽𝐼^∗,𝑎₀₁=−𝛽𝑆^∗,𝑎₁₁=−𝛽,𝑏₁₀=𝛽𝐼^∗,𝑏₀₁=𝑐𝐼^∗/(𝑏+𝐼^∗)²,𝑏₁₁=𝛽,𝑏₀₂=𝑏𝑐/(𝑏+𝐼^∗)³,𝑏₀₃=−𝑏𝑐/(𝑏+𝐼^∗)⁴ and 𝑂(𝑦⁴) 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 𝜎=−3𝜋2𝑎₀₁Δ^3/2×𝑎₁₀𝑏₁₀𝑎²₁₁+𝑎₁₁𝑏₀₂+𝑎₁₀𝑎₀₁𝑏²₁₁+𝑎₁₁𝑏₀₂−2𝑎₁₀𝑏₁₀𝑏²₀₂+𝑎₀₁𝑏₁₀−2𝑎²₁₀𝑏₁₁𝑏₀₂−3𝑏₁₀𝑏₀₃𝑎²₁₀+𝑎₀₁𝑏₁₀, (29) where Δ=𝑎₁₀𝑏₀₁−𝑎₀₁𝑏₁₀=𝐼^∗𝐺(𝐼^∗,𝑐)>0.

By numerical simulation, if we fix 𝑑=0.02,𝜈=26, we know there exist parameter values (𝐴,𝑏,𝑐,𝛽)=(1050,500,128228,0.5), which satisfies ℝ₀=92.93>1 and 𝐻(𝐼^∗,𝑐)=0 such that 𝜎=−97.09<0. On the other hand, there exist values (𝐴,𝑏,𝑐,𝛽)=(1050,0.1,357.3097,0.5), which satisfy ℝ₀=7.29>1 and 𝐻(𝐼^∗,𝑐)=0 too, but 𝜎=0.0015>0. Therefore, there exists an open set 𝑉₁ in the parameter space (𝐴,𝑏,𝑐,𝛽), such that 𝐻(𝐼^∗,𝑐)=0,ℝ₀>1 and 𝜎<0; that is, 𝑉₁=𝐴,𝑏,𝑐,𝛽∶𝐴>0,𝑏>0,𝑐>0,𝛽>0,ℝ₀𝐻𝐼>1,_∗,𝑐=0,𝜎<0. (30) And there exists another open set 𝑉₂ in the parameter space (𝐴,𝑏,𝑐,𝛽), such that 𝐻(𝐼^∗,𝑐)=0,ℝ₀>1 and 𝜎>0; that is, 𝑉₂=𝐴,𝑏,𝑐,𝛽∶𝐴>0,𝑏>0,𝑐>0,𝛽>0,ℝ₀𝐻𝐼>1,_∗,𝑐=0,𝜎>0. (31) By [13], we have the following theorem.

Theorem 6. Suppose ℝ₀>1 and there exists 𝑐_𝑘>0 such that 𝐻(𝐼^∗,𝑐_𝑘)=0 and (𝜕𝐻(𝐼^∗,𝑐)/𝜕𝑐)|_{𝑐=𝑐𝑘}≠0 for the endemic equilibrium 𝐸^∗ of the model. If 𝜎≠0, then a curve of periodic solutions bifurcates from 𝐸^∗ such that the following happens. (1) Suppose (𝐴,𝑏,𝑐_𝑘,𝛽)∈𝑉₁; the model undergoes a supercritical Hopf bifurcation if (𝜕𝐻(𝐼^∗,𝑐)/𝜕𝑐)|_{𝑐=𝑐𝑘}<0 and backward supercritical Hopf bifurcation if (𝜕𝐻(𝐼^∗,𝑐)/𝜕𝑐)|_{𝑐=𝑐𝑘}>0. (2) Suppose (𝐴,𝑏,𝑐_𝑘,𝛽)∈𝑉₂; the model undergoes a subcritical Hopf bifurcation if (𝜕𝐻(𝐼^∗,𝑐)/𝜕𝑐)|_{𝑐=𝑐𝑘}<0 and backward subcritical Hopf bifurcation if (𝜕𝐻(𝐼^∗,𝑐)/𝜕𝑐)|_{𝑐=𝑐𝑘}>0.

4.2. Case 2: ℝ₀<ℝ₀<1

When ℝ₀<ℝ₀<1, model (6) has two endemic equilibria 𝐸_∗(𝑆_∗,𝐼_∗) and 𝐸^∗(𝑆^∗,𝐼^∗) when 𝑎₁<0(𝐴>𝐴₁(𝑐)). Their coordinates satisfy (7) and 𝐼_∗<𝐼^∗. 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 𝐸_∗ and 𝐸^∗ are the same as (24) if one replaces 𝐼_∗ with 𝐼^∗. We have the following result.

Theorem 7. Suppose ℝ₀<ℝ₀<1 and 𝑎₁<0. Then the endemic equilibrium 𝐸_∗ of (6) is a saddle; 𝐸^∗ is a stable node or focus when 𝐻(𝐼^∗,𝑐)>0;𝐸^∗ is an unstable node or focus when 𝐻(𝐼^∗,𝑐)<0;𝐸^∗ is a center of the linear system when 𝐻(𝐼^∗,𝑐)=0.

Proof. Because 𝑑𝑔(𝐼_∗)/𝑑𝐼−𝑑ℎ(𝐼_∗)/𝑑𝐼<0,𝑑𝑔(𝐼^∗)/𝑑𝐼−𝑑ℎ(𝐼^∗)/𝑑𝐼>0, we have 𝐺(𝐼_∗,𝑐)<0,𝐺(𝐼^∗,𝑐)>0 from (27). Then 𝐸_∗ is a saddle, 𝐸^∗ is a stable node or focus when 𝐻(𝐼^∗,𝑐)>0,𝐸^∗ is an unstable node or focus when 𝐻(𝐼^∗,𝑐)<0,and𝐸^∗ is a center of the linear system when 𝐻(𝐼^∗,𝑐)=0.

Here we give some numerical simulations to show that 𝐻(𝐼^∗,𝑐) may be positive or negative for different parameters. For example, we fix 𝐴=15000,𝛽=0.0005,𝑑=0.02,𝑏=0.01,𝜈=34. When 𝑐=7.16, we have ℝ₀=0.50. System (6) has two endemic equilibria, 𝐸_∗(𝑆_∗,𝐼_∗)=(749793.72,0.011)and𝐸^∗(𝑆^∗,𝐼^∗)=(68075.74,400.69). We get 𝐻(𝐼^∗,𝑐)=0.20247. From Figure 3(c), 𝐸^∗ is locally stable. When 𝑐=201, we have ℝ₀=0.019. System (6) has two endemic equilibria too, 𝐸_∗(𝑆_∗,𝐼_∗)=(739115.82,0.589),and𝐸^∗(𝑆^∗,𝐼^∗)=(69059.22,394.41). In this case, 𝐻(𝐼^∗,𝑐)=−0.29239. From Figure 3(d), we can see that 𝐸^∗ is unstable.

From Theorem 7, we know that the endemic equilibrium 𝐸^∗ of system (6) is a center-type nonhyperbolic equilibrium when 𝐻(𝐼^∗,𝑐)=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 𝐸^∗ as we did in Case 1.

Denote nonempty set 𝑉₃ and 𝑉₄ as 𝑉₃=ℝ𝐴,𝑏,𝑐,𝛽∶𝐴>0,𝑏>0,𝑐>0,𝛽>0,₀(𝑐)<ℝ₀<1,𝐻𝐼^∗,𝑉,𝑐=0,𝜎<0₄=ℝ𝐴,𝑏,𝑐,𝛽∶𝐴>0,𝑏>0,𝑐>0,𝛽>0,₀(𝑐)<ℝ₀<1,𝐻𝐼^∗,𝑐=0,𝜎>0. (32)(33)

Theorem 8. Suppose ℝ₀<ℝ₀<1,𝑎₁<0(𝐴>𝐴₁(𝑐)) and there exists 𝑐_𝑘>0 such that 𝐻(𝐼^∗,𝑐_𝑘)=0 and (𝜕𝐻(𝐼^∗,𝑐)/𝜕𝑐)|_{𝑐=𝑐𝑘}≠0 for the endemic equilibrium 𝐸^∗ of model (6). If 𝜎≠0, then a curve of periodic solutions bifurcates from 𝐸₂ such that the following happens. (1) Suppose (𝐴,𝑏,𝑐_𝑘,𝛽)∈𝑉₃; then the model undergoes a supercritical Hopf bifurcation if (𝜕𝐻(𝐼^∗,𝑐)/𝜕𝑐)|_{𝑐=𝑐𝑘}<0 and backward supercritical Hopf bifurcation if (𝜕𝐻(𝐼^∗,𝑐)/𝜕𝑐)|_{𝑐=𝑐𝑘}>0. (2) Suppose (𝐴,𝑏,𝑐_𝑘,𝛽)∈𝑉₄; then the model undergoes a subcritical Hopf bifurcation if (𝜕𝐻(𝐼^∗,𝑐)/𝜕𝑐)|_{𝑐=𝑐𝑘}<0 and backward subcritical Hopf bifurcation if (𝜕𝐻(𝐼^∗,𝑐)/𝜕𝑐)|_{𝑐=𝑐𝑘}>0.

According to the aforementioned results, we know that when 𝑑−𝑏𝛽>0 and 𝑐>𝑐^∗, besides the basic reproduction number ℝ₀, there exists a subthreshold condition for model (3). If ℝ₀≥1, model (3) has a unique endemic equilibrium which is a node or focus. If ℝ₀<ℝ₀<1, model (3) has two endemic equilibria; one is a node or focus and the other is a saddle point. If ℝ₀=ℝ₀, 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. If 𝑑−𝑏𝛽>0 and 𝑐>𝑐^∗, model (3) undergoes a backward bifurcation at ℝ₀=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 ℝ₀, 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 ℝ₀>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 ℝ₀<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 ℝ₀ and it is now necessary to reduce ℝ₀ to a value less than ℝ₀.

𝑐, the maximum treatment per unit of time, is related to the maximum treatment capacity in a city or region. Note that ℝ₀ is a monotone decreasing function of 𝑐. Therefore, under the condition of 𝑑−𝑏𝛽≤0, increasing 𝑐 to a value such that ℝ₀<1 is sufficient to eliminate the disease. While under the condition of 𝑑−𝑏𝛽>0, increasing 𝑐 to a value such that ℝ₀<ℝ₀ is necessary. Decreasing 𝑐 to a value at which ℝ₀<ℝ₀<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 𝐸^∗ or a periodic orbit around 𝐸^∗. 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 𝑐 is small enough such that ℝ₀>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.

(a) The time course of 𝐼(𝑡) when 𝑆(0)=68075,𝐼(0)=400

(b) The time course of 𝐼(𝑡) when 𝑆(0)=68075,𝐼(0)=3

Figure 5 

Different initial states bring different results. The parameters are 𝐴=15000,𝛽=0.0005,𝑑=0.02,𝑏=0.01,𝜈=34,𝑐=7.16, and ℝ0=0.50. The unit of time 𝑡 is day.

The 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 ℝ₀<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 ℝ₀ 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.

Acknowledgments

The authors are grateful to Professor Huaiping Zhu for his valuable comments and suggestions. Hui Wan was supported by NSFC (no. 11201236 and no. 11271196) and the NSF of the Jiangsu Higher Education Committee of China (no. 11KJA110001 and no. 12KJB110012). Jing-an Cui was supported by NSFC (no. 11071011) and Funding Project for Academic Human Resources Development in Institutions of Higher Learning Under the Jurisdiction of Beijing Municipality (no. PHR201107123).