Complex Dynamics of an SIR Epidemic Model with Nonlinear Saturate Incidence and Recovery Rate

Susceptible-infectious-removed (SIR) epidemic models are proposed to consider the impact of available resources of the public health care system in terms of the number of hospital beds. Both the incidence rate and the recovery rate are considered as nonlinear functions of the number of infectious individuals, and the recovery rate incorporates the influence of the number of hospital beds. It is shown that backward bifurcation and saddle-node bifurcation may occur when the number of hospital beds is insufficient. In such cases, it is critical to prepare an appropriate amount of hospital beds because only reducing the basic reproduction number less than unity is not enough to eradicate the disease. When the basic reproduction number is larger than unity, the model may undergo forward bifurcation and Hopf bifurcation. The increasing hospital beds can decrease the infectious individuals. However, it is useless to eliminate the disease. Therefore, maintaining enough hospital beds is important for the prevention and control of the infectious disease. Numerical simulations are presented to illustrate and complement the theoretical analysis.


Introduction
Classical susceptible-infectious-removed (SIR) epidemic models with bilinear incidence rate typically have at most one endemic equilibrium; the disease will die out when the basic reproduction number is less than unity and will persist otherwise [1][2][3][4].However, in practice, many infectious diseases exhibit multiple peaks or periodic oscillations during the outbreak.Therefore, various nonlinear incidence rates have been proposed recently due to the fact that they can produce rich dynamics for the epidemic models.
Liu et al. [5] used the following form of nonlinear saturated incidence rate to incorporate the effect of behavioral changes where βI l S measures the infection force of the disease, 1/(1 + αI h ) describes the inhibition effect from the behavioral change of the susceptible individuals when the number of infectious individuals increases, l, h, β are positive constants and α is non-negative constant [6,7].The nonlinear function g(I) given in (1) has three types, for details one can see [7].Capasso and Serio [8] used the case when l = h = 1, i.e., g(I) = βI 1+αI , to represent a "crowding effect" or "protection measure" in investigating the cholera epidemic in Bari in 1973.Due to the nonlinearity and saturation property of these incidence rates, SIR epidemic models usually possess multiple endemic equilibria and rich nonlinear dynamics [5][6][7][8][9][10][11][12].Furthermore, a compartmental model with nonlinear incidence rate is usually used to explore the impact of intervention strategies on the transmission dynamics of infectious diseases.Therefore, it is essential to investigate the dynamics of this type of epidemic model to prevent and control the spread of infectious diseases.
On the other hand, the resources of health system availability to the public determines how well the diseases are controlled.Particularly, the capacity of the hospital settings and effectiveness and efficiency of the treatment may influence the recovery rate [13].In classical epidemic models, the recovery rate is usually assumed to be proportional to the number of infected, which means that the resources of the health system are quite sufficient for the infectious disease [1].In fact, the number of health care workers and the facilities of the hospital including medical apparatus and equipment, the number of hospital beds and medicines available to the public are limited, especially during the outbreak of the disease.For instance, during the severe acute respiratory syndrome (SARS) outbreaks in 2003, the Chinese government had to create the first and only SARS hospital, Beijing Xiaotangshan Hospital, to treat the larger number of SARS patients, as the normal public-health system and capacity in Beijing City were unable to cope with the rapidly increasing number of SARS cases [9,14].Hence, the capacity of the health care system from both a modeling and analysis of the view should be considered.
Based on the World Health Organization (WHO) Statistical Information System, hospital bed-population ratio (HBPR), the number of available hospital beds per 10,000 population, is used by health planners as a method of estimating resource availability to the public [13,15].To study the impact of HBPR, Shan and Zhu [13] proposed the following nonlinear recovery rate function of the number of hospital beds per 10,000 population b and the number of infective individuals where µ 0 , µ 1 (µ 1 > µ 0 > 0) are, respectively, the minimum and maximum per capita recovery rates.Parameter b is considered as a measure of available hospital resource.Their study showed that the SIR model with standard saturate incidence and recovery rate (2) has rich and interesting dynamics such as backward bifurcation, saddle-node bifurcation, Hopf bifurcation and cusp type of Bogdanov-Takens bifurcation of codimension 3. Recently, Abdelrazec et al. [16] applied the recovery rate (2) to investigate the impact of available resources of the health system on the spread and control of dengue fever, which could be helpful for public health authorities in their planning of a proper resource allocation for the control of dengue transmission.Motivated by these points, our model thus incorporates both nonlinear incidence rate and recovery rate to well control the emerging infectious.In other words, the combined effects of government intervention and hospitalization condition are considered to prevent an outbreak.However, it is more natural to consider perturbations of contagion coefficients through the Wiener process or treat them directly as random variables since the transmission coefficients are usually unknown in practice (see, for example, [17,18] and references therein).Here, we focus on the deterministic epidemic model and leave the consideration of randomness in epidemiological model as our future work.Therefore, in this paper, we investigate the following deterministic epidemic model with a nonlinear incidence rate and recovery rate: with µ(b, I) defined in (2).In system (3), S(t), I(t) and R(t) are the numbers of susceptible, infectious and recovered individuals at time t, respectively; A > 0 is the recruitment rate of the population; d > 0 is the per capital natural death rate of the population; ≥ 0 is the per capita disease-induced death rate; µ(b, I) is the per capita recovery rate of infectious individuals incorporating the impact of the capacity and limited resources of the health care system, and b > 0 is the number of available hospital beds per 10,000 population.The organization of this paper is as follows.In the next section, we investigate the existence and classification of the equilibria for system (3).In Section 3, we analyze the nonlinear dynamics for system (3) such as forward bifurcation, backward bifurcation, saddle-node bifurcation and Hopf bifurcation.In Section 4, the numerical simulations are obtained to verify our results.A brief discussion is then presented and conclusions are presented in the final section.

Existence and Classification of Equilibria
Since the first two equations in system (3) are independent of variable R, it suffices to consider the following reduced system: It should be noted that the total population number N = S + I + R satisfies the equation Therefore, the biologically feasible region is a positively invariant with respect to model (4).

Existence of Equilibria
It is obviously that system (4) always admits a disease-free equilibrium E 0 (A/d, 0).Following the techniques of van den Driessche and Watmough [19], the basic reproduction number of system (4) can be expressed as The endemic equilibrium of system (4) can be obtained by solving the following algebraic equations From the second equation of ( 6), we have Substituting (7) into the first equation of ( 6) leads to where where δ 0 = d + + µ 0 and δ 1 = d + + µ 1 .We will use these two notations in the whole paper.If ∆ 0 > 0, the roots of (8) read If ∆ 0 = 0, the two roots I 1 and I 2 of (8) coalesce into a unique root with multiplicity 2 denoted as I * = − a 1 2a 0 .If I i > 0, system owns endemic equilibrium E i (S i , I i ) with If I * > 0, the two endemic equilibria E 1 and E 2 coalesce into one endemic equilibrium E * (S * , I * ), where Based on the above statements, we investigate the existence of equilibria in the following three cases.
In this case, we have a 0 > 0, a 2 = 0 and ∆ 0 = a 2 1 > 0. Quadratic Equation (8) has no positive root if a 1 ≥ 0 and a unique positive root I 2 if a 1 < 0. One can verify that a 1 < 0 is equivalent to In the case of guarantees the inequality (10) holds.Then, system (4) has a unique endemic equilibrium E 2 .
In this case, we have that a 0 > 0 and a 2 > 0. It is easy to show that (8) has no positive root if a 1 ≥ 0, and inequality (10) implies that , then a 1 < 0. The number of roots for Equation (8) determined by ∆ 0 , namely, Equation ( 8) has no positive root if ∆ 0 < 0, one root I * of multiplicity 2 if ∆ 0 = 0, and two positive roots I 1 and I 2 if ∆ 0 > 0. Next, we determine the signs for ∆ 0 .By considering ∆ 0 as a quadratic equation of variable 1 − R 0 , a straightforward computation derives that ∆ 0 > 0 when R 0 > P 1 or R 0 < P 2 ; ∆ 0 = 0 when R 0 = P i , (i = 1, 2), and ∆ 0 < 0 when P 2 < R 0 < P 1 , where and Notice that 0 < P 2 < P 0 < P 1 < 1.Therefore, system (4) has no endemic equilibrium if R 0 < P 1 , one endemic equilibrium E * of multiplicity 2 if R 0 = P 1 , and two endemic equilibria E 1 and E 2 if Summarizing the discussions above, we have the following results.

1.
The disease-free equilibrium E 0 always exists.
If R 0 < 1, and , there is no endemic equilibrium; , system (4) has two endemic equilibria E 1 and E 2 when P 1 < R 0 < 1, and these two equilibria coalesce into one endemic equilibrium E * when R 0 = P 1 ; otherwise, there is no endemic equilibrium.

Types of Equilibria
In this section, we investigate the stabilities of the equilibria for system (4).Let E(S, I) be any equilibrium of system (4).Then, the Jacobian matrix J(E) of system (4) around E(S, I) can be expressed as The corresponding characteristic equation of J(E) is where TrJ(E) and DetJ(E) are the trace and the determinant of matrix J(E), respectively.
Theorem 2. For system (4), Proof.Direct calculation yields that −d and (d + + µ 1 )(R 0 − 1) are the two roots of the Jacobian matrix J(E 0 ) at disease-free equilibrium E 0 .Then, E 0 is an attracting node if R 0 < 1, while it is a hyperbolic saddle if R 0 > 1.
If R 0 = 1, the second eigenvalue is zero.To determine the type of E 0 , we first transform the disease-free equilibrium E 0 of system (4) into the origin.Let U = I and V = S − A/d; then, we have Expanding system (13) in Taylor series at (U, V) = (0, 0) to the second order, it follows that Using the transformation then we can rewrite system (14) as , it is unnecessary to calculate the center manifold, and system (15) already shows that E 0 is a saddle-node.
if δ = 0, E 2 is a weak focus or a center; 3.
if δ > 0, E 2 is a attracting node or focus.
Proof.Denoting the endemic equilibrium point as Ē( S, Ī), the corresponding characteristic equation of J( Ē) is where and We prove the sign of DetJ( Ē) in two cases.
Rewriting (18) as follows: where Since I 2 satisfies Equation ( 8), then we can simplify ϕ(I 2 ) as where Using Equations ( 21) and ( 22), we know that TrJ(E 2 ) and δ have the opposite sign.Thus, if δ < 0, E 2 is a repelling node or focus; if δ > 0, E 2 is an attracting node or focus; and if δ = 0, E 2 could be a weak focus or a center.
, then the unique disease-free equilibrium E 0 is globally asymptotically stable, whereas if , then the endemic equilibrium E 2 is globally asymptotically stable.
Proof.If R 0 = P 1 , it follows from Theorem 1 (iv) that two endemic equilibria E 1 and E 2 coalesce into a unique E * .One can easily obtain that 0 and TrJ(E * ) are the two eigenvalues of characteristic equation for J(E * ).Next, we only consider the case TrJ(E * ) = 0. Similar to the computational process in Theorem 2, we linearize system (4) at E * and diagonalize the linear part.After complication calculation, we can transform system (4) as the following system: Obviously, E * is a saddle-node point.Therefore, according to Chapter 20.1C in [20] and Theorem 1 (iv), we know that system (4) undergoes a saddle-node bifurcation when R 0 passes through the critical value P 1 .This proof is completed.

Hopf Bifurcation
Next, we study the Hopf bifurcation of system (4).From the above discussion, Hopf bifurcation can only occur at endemic equilibrium E 2 and the necessary condition of Hopf bifurcation requires that TrJ(E 2 ) = 0.The proof of Theorem 3 implies that TrJ(E 2 ) = 0 if and only if δ = 0, and DetJ(E 2 ) > 0 if E 2 exists.Thus, the eigenvalues of J(E 2 ) have a pair of pure imaginary roots if and only if It follows from Theorem 3.4.2 in [21] that system (4) maybe undergo a Hopf bifurcation at endemic equilibrium E 2 if δ = 0.
Next, we study the normal form of system (4).Let u = S − S 2 and v = I − I 2 ; then, system (4) becomes Expanding system (29) in Taylor series at (u, v) = (0, 0) to the third order, we have where and Using the translation we can rewrite system (30) as Notice that system (31) is exactly the normal form of system (4).Based on Theorem 3.4.2 in [21], the discriminatory quantity σ is given by where f xy denotes (∂ 2 f /∂x∂y)(0, 0), etc.Thus, after an extensive calculation, we have σ = 1 8a 3  12 ω 2 − a 12 (a If the discriminatory quantity σ is not zero, we have the following result. Theorem 7. System (4) undergoes a Hopf bifurcation at endemic equilibrium E 2 if δ = 0.Moreover, if σ < 0 (resp., σ > 0), then at least one attracting (resp., repelling) limit cycle bifurcates from E 2 .

Bifurcation Diagram and Simulation
The theoretical results in previous sections show that the dynamics of system (4) depend on the expressions of R 0 , a 1 , ∆ 0 and δ.Therefore, in the following, we first analyze these expressions in the From the expression of the basic reproduction number R 0 , we know that R 0 = 1 defines a straight line L 0 in µ 1 − b plane One can easily verify that hyperbola C a and line ), and f a (µ 1 ) is a decreasing convex function of µ 1 with a horizontal asymptote b = 0.
If ∆ 0 (µ 1 , b) = 0, we obtain the curves C ± ∆ by solving b in term of µ 1 , where Therefore, the curve C − ∆ under the curve C a is convex and decreasing with the asymptote b = 0.
where m 1 , m 2 are given in (23).Furthermore, we define which will be used later.Based on the above discussions and Theorem 1, let Then, system (4) has no endemic equilibrium in region Ω 0 , one endemic equilibrium E 2 in region Ω 1 and two equilibria E 1 and E 2 in region Ω 2 .The two equilibria in Ω 2 coalesce into one equilibrium E * on curve C − ∆ .For more intuitive observation, one can see Figure 1a.The stability of these equilibria can be observed in Figure 1b.According to Theorem 2, the disease-free equilibrium E 0 always exists, and it is locally asymptotically stable (L.A.S) in region Ω 0 ∪ Ω 2 and unstable in region Ω 1 .Furthermore, Theorem 4 suggests that E 0 is globally asymptotically stable (G.A.S) in region I (green region in Figure 1b).From Theorem 3, we know that E 1 is always unstable whenever it exists, and E 2 is unstable in region II (red region in Figure 1b) and is L.A.S. in region II 1 (light blue region in Figure 1b).Moreover, Theorem 4 implies that E 2 is G.A.S in region II 2 (blue region in Figure 1b).Finally, we choose a point (marked in black dot in Figure 1b) in the region III and plot its phase portrait in Figure 2. It is clear that the disease-free equilibrium E 0 and two endemic equilibria E 1 and E 2 coexist, and E 0 is L.A.S, E 1 is a saddle and E 2 is unstable.According to Theorem 5, system (4) undergoes forward bifurcation on L + 0 , backward bifurcation on L − 0 , pitchfork bifurcation when transversally passing through the line L 0 at point b * and saddle-node (SN) bifurcation on the curve C − ∆ when the two endemic equilibria coalesce into one endemic equilibrium E * .Hopf bifurcation occurs on the curve C δ as proved in Theorem 7. Figure 3a shows the Hopf bifurcation curve in the µ 1 − b plane.We choose one point (µ 1 , b) = (36, 0.615) marked with a black asterisk below the blue curve in Figure 3a, and plot the phase portrait at this point in Figure 3b.We can observe from Figure 3b that the trajectory (blue curve) starting at (184.2, 0.042) spirals outward to the stable limit curve (black curve) and the trajectory (red curve) starting at (188, 0.1) spirals inward to the stable limit curve.At point (µ 1 , b) = (36, 0.615), one can easily obtain that E 1 does not exist and E 2 exists, but it is unstable.Therefore, system (4) has a stable limit curve that enriches the unstable equilibrium E 2 .The typical bifurcation diagrams in R 0 − I plane can be seen in Figure 4. Figure 4a shows that there is a forward bifurcation at R 0 = 1 from disease-free equilibrium E 0 to a unique endemic equilibrium E 2 when R 0 > 1.If we decrease b from 0.65 to 0.61, system (4) not only undergoes forward bifurcation but also undergoes Hopf bifurcation as shown in Figure 4b.In this case, a stable limit cycle bifurcated from forward Hopf bifurcation and disappears from the backward Hopf bifurcation.If we further decrease b to 0.07, we can observe from Figure 4c that the backward bifurcation and saddle-node bifurcation occur.This illustrates that the number of hospital beds plays an important role in controlling an infectious disease.Finally, we investigate the effect of the nonlinear incidence rate.For the nonlinear saturate incidence rate βSI/(1 + αI), it is known that larger α reflects stronger inhibition effect caused by the infective individuals.Although the increasing of α cannot reduce the basic reproduction number R 0 , it does affect the dynamical behaviors of the system.Based on previous results, we can see from Figure 6a that backward bifurcation occurs on green lines, saddle-node bifurcation occurs on red curves and Hopf bifurcation happens on blue curves.Figure 6a and the magnified diagram Figure 6b also illustrate that the occurrence of backward bifurcation, saddle-node bifurcation and Hopf bifurcation shrink with increasing α.That is to say, the strong inhibition effect will lower the complexity of the transmission of the infectious diseases.In other words, enhancing the public's defensive "crowding effect" through education or medium is beneficial to control and prevent the spread of the infectious diseases.Furthermore, Figure 6 also implies that the stronger inhibition effect of the lower number of hospital beds will be needed to eliminate the diseases.It reveals that strengthened public defensive "crowding effect" can mitigate the demand for hospital beds in controlling the spread of the disease during an outbreak.

Discussion
Due to the important biological significance of hospital beds, in this paper, we have investigated an SIR epidemic model to simulate the impact of a limited health care system in terms of the number of hospital beds.Theoretical analysis and numerical simulations illustrated that system (4) has at most three equilibria and possesses complex nonlinear dynamics.From Theorems 5-7 and Figures 4 and 5, we have shown that system (4) can undergo backward bifurcation, saddle-node bifurcation and Hopf bifurcation due to the insufficient number of hospital beds.
If R 0 > 1, it follows from Theorem 3 and Figure 1b that E 2 is the unique endemic equilibrium, and it is stable if b is sufficiently large such that the values of b are all above the curve C δ .Since I 2 > 0 and ∂I 2 (b) ∂b | b>0 < 0, then I 2 (b) is a monotone decreasing function of b and I 2 (b) trends to some positive constant as b to infinity.That is, increasing the number of hospital beds can only reduce the number of the infectious, but cannot eliminate the infectious disease (see Figure 5a,b).If R 0 < 1, it follows from Theorem 3 and Figure 1b that the disease can be eliminated when the value of b above the curve is C − ∆ .This means that the basic reproduction number is not the only evaluation standard for the control and elimination of the disease.It also depends on the resources of the health care system such as the number of hospital beds.
The bifurcation analysis is carried out through reducing the three dimension model (3) to the planar system (4).The typical bifurcation curves and diagrams are shown in Section 4, and all of the results reveal that system (3) possesses complex dynamics such as forward bifurcation, backward bifurcation, saddle-node bifurcation and Hopf bifurcation.Contrasting to the results for the classic SIR epidemic models, we also find that the nonlinearity of incidence rate and recovery rate are important factors that lead to very rich dynamics.Moreover, we find that the dynamical behavior of system (3) not only depends on R 0 but also relies on the number of hospital beds and the inhibition effect caused by the infective individuals.Therefore, the public health makers should consider the combined effects of the government intervention strategies (such as education and medium) and hospitalization conditions to make guidelines in controlling the spread of infectious diseases.Hopefully, in the future, we can explore more theoretical results to control and eliminate infectious diseases, especially in the outbreak of diseases.
use the Dulac criteria to prove the global stability of E 2 if R 0 > 1. Define the Dulac function as B(S, I) = 1 + αI I , and let f = A − βSI 1 + αI − dS and g = βSI 1 + αI − dI − I − µ(b, I)I.

Figure 1 .
Figure 1.(a) the bifurcation curve in the (µ 1 , b) plane; and (b) the stability curve of equilibria in the (µ 1 , b) plane.

Figure 3 .
Figure 3. (a) Hopf bifurcation curve in µ 1 − b plane.The black asterisk below the curve is the point we choose to plot the portrait in (b).Here, other parameters have the same values as in Figure 2.

Figure 4 .
Figure 4. Bifurcation diagrams of system (4) in R 0 − I plane with different b.The blue curves represent the stable fix points or limit cycles and the red dashed curves represent the unstable fix points.To further explore the impact of the number of hospital beds, we also present the possible bifurcation diagrams in b − I plane by fixing all the parameters except b.Figure 5 depicts some typical bifurcation diagrams in b − I plane with different R 0 .If R 0 > 1, one can see from Figure 5a,b that increasing the beds can reduce the number of infectious individuals but can not eliminate the disease.Especially in the case of R 0 = 1.29,Hopf bifurcation occurs and a stable limit cycle bifurcated from forward Hopf bifurcation disappears from backward bifurcation.If R 0 < 1, backward bifurcation and saddle-node bifurcation may occur and the disease dies out when the value of b is above the curve C −

Figure 5
depicts some typical bifurcation diagrams in b − I plane with different R 0 .If R 0 > 1, one can see from Figure 5a,b that increasing the beds can reduce the number of infectious individuals but can not eliminate the disease.Especially in the case of R 0 = 1.29,Hopf bifurcation occurs and a stable limit cycle bifurcated from forward Hopf bifurcation disappears from backward bifurcation.If R 0 < 1, backward bifurcation and saddle-node bifurcation may occur and the disease dies out when the value of b is above the curve C − ∆ based on Theorem 4 and Figure 1.

Figure 5 .
Figure 5. Bifurcation diagrams of system (4) in b − I with different R 0 .

Figure 6 .
Figure 6.(a) impact of the nonlinear incidence rate on dynamical behaviors of system (3); (b) the magnified diagram of black box marked in figure (a).Besides α, other parameters have the same values as in Figure 2.