Dynamic behaviors of a modified SIR model in epidemic diseases using nonlinear incidence and recovery rates

The transmission of infectious diseases has been studied by mathematical methods since 1760s, among which SIR model shows its advantage in its epidemiological description of spread mechanisms. Here we established a modified SIR model with nonlinear incidence and recovery rates, to understand the influence by any government intervention and hospitalization condition variation in the spread of diseases. By analyzing the existence and stability of the equilibria, we found that the basic reproduction number R0 is not a threshold parameter, and our model undergoes backward bifurcation when there is limited number of hospital beds. When the saturated coefficient a is set to zero, it is discovered that the model undergoes the Saddle-Node bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation of codimension 2. The bifurcation diagram can further be drawn near the cusp type of the Bogdanov-Takens bifurcation of codimension 3 by numerical simulation. We also found a critical value of the hospital beds bc at R0<1 and sufficiently small a, which suggests that the disease can be eliminated at the hospitals where the number of beds is larger than bc. The same dynamic behaviors exist even when a ≠ 0. Therefore, it can be concluded that a sufficient number of the beds is critical to control the epidemic.


Introduction
Since the development of the first dynamic model of smallpox by Bernoulli in 1760, various mathematical models have been employed to study infectious diseases [1] in order to reveal the underlying spread mechanisms that influence the transmission and control of these diseases. Among them, Kermack and Mckendrick [2] initiated a famous SIR type of compartmental model in 1927 for the plague studies in Mumbai, and succeed in unveiling its epidemiological transition. Since then, mathematical modeling has become an important tool to study the transmission and spread of epidemic diseases.
In the modeling of infectious diseases, the incidence function is one of the important factors to decide the dynamics of epidemic models. Bilinear and standard incidence rates, both monotonically increasing functions of the total of infected class, have been frequently used in early PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 1 / 28 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 epidemic models [3]. In those models, the dynamics of models are relatively simple and almost determined by the basic reproduction number R 0 : the disease will be eliminated if R 0 < 1, otherwise, the disease will persist. However, intervention strategies, such as isolation, quarantine, mask-wearing and medical report about emerging infectious diseases, play an vital role in controlling the spread, sometimes contributing to the eradication of diseases. For instance, the SARS in 2003 and novel influenza pandemic in 2009 have been well controlled by taking these intervention actions [4][5][6][7][8][9][10][11][12][13][14][15]. Hence, it is essential to expand the modeling studies to the investigation of the combined effects of these major intervention strategies. The generalized models will provide further understanding of the transmission mechanisms, and modify guidelines for public health in control of the spread of infectious diseases.
In recent years, a number of compartmental models have been formulated to explore the impact of intervention strategies on the transmission dynamics of infectious diseases. If denote the total number of hospital individuals, exposed and infectious as N, E and I respectively, Liu et al [16], Cui et al [17,18] used βe −mI , be À a 1 EÀ a 2 IÀ a 3 H and c 1 − c 2 f(I) to study the impact of media coverage on the dynamics of epidemic models, respectively. However, people may adjust their behaviors according to these government intervention. Therefore, Ruan and Xiao [19] set incidence function in the form of f 1 ðIÞS ¼ KIS 1þaI 2 (a special case of kI p S 1þaI q [20,21]) to include the above "psychological" effect: when the number of infectious individuals increases and is reported through social media, the susceptible individuals will stay alert spontaneously to reduce any unnecessary contact with others, thus lowering the contact and transmission incidence.
On the other hand, medical treatments, determining how well the diseases are controlled, are normally expressed as constant recovery rates in the current models. These recovery rates depend on various health systems and hospitalization conditions, such as the capacity of the hospitals and effectiveness of the treatments. Advanced models (see [22][23][24]) started to corporate the limited medical resources into the spread dynamics of infectious diseases. In the literature [22], Wang and Ruan first introduced a piece-wise treatment function in an SIR model, where the maximal treatment capacity was used to cure infectives so that the epidemic of disease can be controlled. This situation occurs only if the infectious disease needs to be eliminated due to its threats to public. They discovered that the model undergoes Saddle-Node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation, standing for the collision of two equilibria, the existance of periodic diseases, and two varying parameters in system, respectively. Wang [23] further modified the treatment rate to be proportinal to the number of infectives before the capacity of hospital was reached, by The model was then found to perform backward bifurcation [23], indicating that the basic reproduction number was no longer a threshold.
In common hospital settings, the number of beds is an indicator of health resources, particularly the medical treatments of the infectives. Under this consideration, Shan and Zhu [24] defined the recovery rate as a function of b, the number of hospital beds, and I, the number of infectives.
where μ 0 is the minimum per capita recovery rate, and μ 1 the maximum per capita recovery rate. They chose the standard incidence rate and discovered the complicated dynamics including Saddle-Node bifurcation, Backward bifurcation and Bogdanov-Takens bifurcation of codimension 3, which means that the recovery rate contributes to the rich dynamics of epidemic models. Our strategy thus becomes, both government intervention and hospitalization condition need to be incorporated to achieve a better control of the emerging infectious. Therefore, the incidence rate is expressed as where a is positive constant and c > À 2 ffiffi ffi a p (so that aI 2 + cI + 1 > 0 for all I > 0 and hence β(I) > 0 for all I > 0). When the threshold of the number of infected individuals I Ã is reached, the contact transmission rate starts to decrease as the number of infected individuals grows. As shown in Fig 1, the incidence β(I) increases to its maximum and then decreases to zero as I tends to infinity, which explains the phenomenon where the rate of contacting between infected I and susceptible S decreases after government intervention. We use the same expression of hospitalization conditions as the literature [24], and the following model is then established, where A is the recruitment rate of the susceptible population, d the natural death of the population, and α the per capita disease-induced death rate, respectively. For system (3), the cone R 3+ is a positive invariant. The C 1 smoothness of the right side of system (3) implies the local existence and uniqueness of the solution with initial values in R 3+ . Adding up the three equations of system (3), we get N 0 (t) = A − dN − αI. Therefore, all solutions in the first octant approach, entering or staying inside the set, are defined by This paper will be organized as follows. In section 2, we study the existence of the equilibria of our model. In section 3, we study the stability of the equilibria. In section 4, we examine the dynamics of the model by first looking at the backward bifurcation of system, then the much complicated Hopf bifurcation and Bogdanov-Takens bifurcation of codimension 2 and 3. We summarize our results and discuss the epidemiological significance of the number of hospital beds and intervention strategies in section 5.

Existence of equilibria
For simplicity we will focus on the case when c = 0. If c 6 ¼ 0, but c is in small neighborhood of zero, the behaviors of model still exist. Our model thus becomes, Since the first two equations are independent of the third, it suffices to consider the first two equations. Thus, we will focus on the reduced model in the following discussions.
We find equilibria by setting the right hand of system (5) equal to zero: Obviously, a trivial solution of Eq (6) is E 0 ðS; IÞ ¼ A d ; 0 À Á , a disease free equilibrium(DFE). For E 0 , by using the formula in [25], one can calculate the reproduction number For any positive equilibrium E(S, I), also called endemic equilibrium, when exists, its S and I coordinates satisfy where the I coordinate will be the positive root of the following cubic equation Denote Δ 0 the discriminant of f(I) with respect to I, then Case 2: R 0 < 1. In this case, D > 0. If C > 0, equation f(I) = 0 has no positive solution (see Fig 1(a)). If C < 0, similar to Lemma 2.1 described by Huang and Ruan [26], the following conclusions can be drawn as shown by Fig 1(b). (c) Δ 0 > 0, we find that Eq (9) has no positive solution.
(b) C < 0, Δ 0 < 0, there exist two endemic equilibria E 1 ð " I 1 ; " S 1 Þ and E 2 ð " I 2 ; " S 2 Þ, and when Δ 0 = 0 these two endemic equilibria coalesce into the same endemic equilibrium E Ã ; otherwise system (5) has no endemic equilibrium. Remark 0.2. By calculation, we get that When a is sufficiently small and tends to zero, the sign of Δ 0 will be determined by the zero power of a. Therefore,

Stability analysis of equilibria
In order to discuss the stability of equilibrium, we need the Jacobian matrix of system (5) at any equilibrium E(S, I). If we denote the Jacobian as J(E) = (j ij ) 2×2 , i, j = 1, 2, then a straightforward calculation gives Firstly, we present a theorem about the disease-free equilibrium E 0 (A/d, 0). Theorem 0.3. E 0 is an attracting node if R 0 < 1, and hyperbolic saddle if R 0 > 1. When and attracting semi-hyperbolic node For system (5), −d and d 1 ðR 0 À 1Þ are two eigenvalues of J(E 0 ). So, E 0 is an attracting node if R 0 < 1, and unstable if R 0 > 1. When R 0 ¼ 1, the second eigenvalue becomes zero. In order to analyze the behavior of E 0 , we linearize system (5) and use the transformation of where P(X, Y) and Q(X, Y) represent the higher order terms. Obviously, E 0 is a saddle-node if , applying the center manifold theorem, system (5) becomes where Q 1 (X, Y) represents the higher order term. Thus, E 0 is an attracting semi-hyperbolic node of codimension 2.
Theorem 0.4. If dδ 0 > βA, E 0 is globally asymptotically stable. Proof. If dδ 0 > βA, it is obvious that R 0 < 1 and C > 0. From Theorem 0.1 and 0.3, E 0 is the unique attracting node of system (5). In order to prove that the disease free equilibrium E 0 is globally and asymptotically stable, we construct the following Liapunov function: It is easy to discover that E 0 A d ; 0 À Á attains the global minimum of the function V(S, I), so V(S, I) > 0. Along system (5), it turns out: Since The equality _ V ðS; Let Eð " S; " I Þ be any endemic equilibrium, one can verify that its characteristic equation can be written as Obviously, the signs of the eigenvalues are determined by f 0 ð " I Þ and tr(J E ). From Fig 2, we know that f 0 ð " I 1 Þ < 0; f 0 ð " I 2 Þ > 0, so E 1 is a hyperbolic saddle and E 2 is an anti-saddle. E 2 is an attracting node or focus, if tr(J E ) < 0; E 2 is a weak focus or a center, if tr(J E ) = 0; E 2 is a repelling node or focus, if tr(J E ) > 0. So we obtain the following theorem. Theorem 0.5. For system (5), there are two endemic equilibria E 1 , E 2 when R 0 < 1 and Δ 0 < 0. Then the low endemic equilibrium E 1 is a hyperbolic saddle, and the higher endemic equilibrium E 2 is an anti-saddle. When R 0 > 1 there is a unique endemic equilibrium, which is an anti-saddle.
Proof. For convenience of the proof, we suppose that the total number of the population is N(t). System (4) becomes the following system Á T and we can write Eq (20) in vector forms as: where Then, We know that the dominant eigenvalue of H(V 0 ) is zero, if R 0 ¼ 1. It is well known that we can decompose a neighborhood of the disease-free state into stable manifold W s and a center manifold W c . Thus, the dynamic behavior of system (20) can be determined by the flow on the center manifold. We know that zero is a simple eigenvalue and the W c is tangential to the eigenvector V 0 at V 0 . Thus, we can assume that W c has the following form: In other words, W c can be decomposed into two components. The α gives the component of the center manifold that lies along the dominant eigenvector; the component that does not lay along the dominant eigenvector can be given by Z(α). So, V Ã Á Z(α) = 0. In order to determine the dynamic behavior of system (20), we just need to see how α depends on time t. Let since W c is an invariant, from Eq (21) we have Multiplying both sides of the above equation by V Ã and using the following equations: and The sign of this expression for small α is what determines whether the disease can invade at the bifurcation point. In the limit, as α goes to zero, Eq (28) becomes: where which gives the rate of change of the vector field as the disease invades. Hence, the number determines whether the disease can invade when R 0 ¼ 1, and hence gives the sign of the bifurcation. For our system, by computation we can get the V Ã and V 0 as follows: where Then we can get that According to [27], we know that system (5) undergoes backward bifurcation, when h > 0, i.e., b < dðm 1 À m 0 Þ bd1 ; and system (5) undergoes forward bifurcation, when h < 0, i.e., b > dðm 1 À m 0 Þ bd1 . Proposition 0.7. When R 0 passes through R c 0 and tr(I Ã ) 6 ¼ 0, system (5) with a ! 0 undergoes a saddle-node bifurcation.
If tr(I Ã ) 6 ¼ 0, we can linearize system (5) at the E Ã and diagonalize the linear part. Then we can get the following form Where T is the non-singular transformation to diagonalize the linear part.
saddle-node if tr(I Ã ) 6 ¼ 0. Combined with Theorem 0.1, system (5) undergoes saddle-node bifurcation when R 0 passes through the critical value R c 0 , as a ! 0. If tr(I Ã ) = 0, E Ã is a cusp and we will prove it in the next section.
Based on the above analysis, we know that system (5) undergoes some bifurcation. In order to consider the impact of hospital bed number and the incidence rate on the dynamics of the model, we will chose b and β as bifurcation parameters to describe these bifurcations. The basic production number R 0 ¼ 1 defines a straight line C 0 in the (β, b) plane, C ¼ 0 also defines one branch of the hyperbolic C B (see Fig 3), The branch of C B intersects with C 0 at the point K dd and with β-axis at the It is easily found that f C is an increasing convex function of β in the first quadrant. Let g; Here C þ 0 and C À 0 are two branches of C 0 joint at point K.

Dynamic behaviors of a modified SIR model
Define the curve Δ 0 (β, b) = 0 as C D 0 , one can verify that Hence, the curve C D 0 is tangent to the curve C 0 at the point K and the β-axis at the point K 0 , which means that K is the only point at which the curve C D 0 is tangent to the curve C 0 .
If b = 0, D 0 ðb; 0Þ ¼ D 0 ðbÞ ¼ À 3d 0 ðdd 0 À bAÞ 2 g 1 ðbÞ; where Through computing, we find that equation Δ 0 (β, 0) = 0 has three real solutions It is easy to verify β 0 > β 1 , so the curve C D 0 will not intersect with the abscissa axis when Obviously Based on the above the discussion and Theorem 0.1, if we define ; D 0 ðb; bÞ < 0g; then there is one endemic equilibria in the region D 1 and two endemic equilibria in the region D 2 . System (5) undergoes saddle-node bifurcation on the cure C D 0 when b 2 dd 0 A ; dd 1 A À Ã . The backward bifurcation occurs on the C À 0 and forward bifurcation occurs on the C þ 0 . The pitchfork bifurcation occurs when transversally passing through the curve C 0 at the point K. Especially, if a = 0, system (5) has a semi-hyperbolic node of codimension 2 at the point K and we can solve b in term of β from Δ 0 (β, b) = 0, Now, we discuss the Hopf bifurcation of system (5). It follows from Eq (18) that if Hopf bifurcation occurs at one endimic equilibrium Eð " S; " I Þ, we have tr(J E ) = 0. Note that from Eq (19) we can rewrite tr(J E ) as where This means sufficient number of hospital beds excludes the possibility of the disease oscillation. From the expression of h 1 ð " I Þ in Eq (37), it is not an easy task to study the Hopf bifurcation from the polynomial Eq (37), we will study a simple case when a = 0, and give the simulations to explore the case when a is small. When a = 0, the polynomial Eq (37) is reduced to One can verify the following lemma In order to study Hopf bifurcation and Bogdanov-Takens bifurcation, we will assume that b < m 1 À m 0 À 2d b . Denote the discrimination of h 1 ð " I Þ as Δ 1 . Since b 1 < 0, function h 1 ð " I Þ must have one negative real root. As shown in Fig 4, It is not difficult to verify function h 1 ð " IÞ has two humps which locate on the different sides of vertical axis, and the maximum is obtained on the left hump, while the minimum is obtained on the right hump.
The number of roots of function h 1 ð " I Þ is determined by the sign of the Δ 1 . When exist, we denote the roots as H m and H M with H M ! H m . Lemma 0.9. @I 2 @b > 0; 8b > 0 and @I 2 @b < 0; 8b > 0. Proof. From the Eq (9) and the expression of I 2 , direct calculation leads to One can find that @C @b ¼ bd 1 > 0. We will prove @I 2 @b < 0 in two cases. If R 0 < 1, then @D @b ¼ dd 1 À bA > 0. Recall the analysis of the existence of the equilibria we know that the C < 0, so the @I 2 @b < 0. If R 0 > 1, then lim b!þ1 so the @I 2 ðbÞ @b is an increasing function of b with supremum 0 in the (0, +1), so for 8b > 0  From Lemma 0.9, we have @I 2 @b j b¼b < 0, so dg dm 1 < 0 (one can verify that H m < b). Then the proof of theorem is completed.
The reason why we choose different parameters to unfold Hopf bifurcation in Theorem 0.10 is that the transversality condition may fail at some point if we focus on one parameter.
In order to verify that Hopf bifurcation occurs in the system, we need to know the type of E 2 . If E 2 is a weak focus, Hopf bifurcation can happen, otherwise system does not undergo Hopf bifurcation. Because system (5) is analytic when a = 0, E 2 can only be weak focus or center. We can distinguish these two types of singularities by calculating Lyapunov coefficients and verifying it through numerical simulation.
We choose a point(β, b) = (0.3683, 0.1587) in Fig 5(a) to plot the phase portrait at the point. In Fig 5(b), as t ! +1, the trajectory starting at (9, 0.1) spirals outward to the stable limit cycle (red curve) and E 2 (8.0794, 0.19363) is stable. Because system (5) is a plane system, there must exist a unstable limit cycle between the stable endemic equilibria and stable limit cycle (black curve). The blue curve in the Fig 5(b) is the unstable manifold of E 1 .

Bogdanov-Takens bifurcation
From Theorem 0.1 we know that the two equilibria E 1 and E 2 coalesce at the equilibria We can find that det(I Ã ) = 0 in Eq (18) if R 0 ¼ R c 0 . From Proposition 0.7, we know that E Ã is a saddle-node point if tr(I Ã ) 6 ¼ 0. If tr(I Ã ) = 0, Eq (18) has a zero eigenvalue with multiplicity 2, which suggests that system (5) may admit a Bogdanov-Takens bifurcation. Then, we give the following theorem.
Theorem 0.11. For system (5) with a = 0, suppose that Δ = 0, h(I Ã ) = 0 and h 0 (I Ã ) 6 ¼ 0, then E Ã is a Bogdanov-Takens point of codimension 2, and system (5) localized at E Ã is topologically equivalent to Proof. Changing the variables as x = S − S Ã , y = I − I Ã , then system (5) becomes dx dt ¼ À ðd þ bI Ã Þx À bS Ã y À bxy; where tr(I Ã ) = 0 and det(I Ã ) = 0 imply that so the generalized eigenvectors corresponding to λ = 0 of Jacobian matrix in system (5) at the point E Ã are Let T = (T ij ) 2×2 = (V 1 , V 2 ), then under the non-singular linear transformation here Using the near-identity transformation and rewrite u, v into X, Y, we obtain To consider the sign of M 21 , note that For system (5), the condition of the existence of endemic equilibrium is b < dðm 1 À m 0 Þ bd 1 , hence, M 21 < 0. Then we will determine the sign of M 22 by If h 0 (I Ã ) 6 ¼ 0, we make a change of coordinates and time and preserve the orientation by time In [28], a generic unfolding with the parameters ε = (ε 1 , ε 2 , ε 3 ) of the codimension 3 cusp singularity is C 1 equivalent to About system (47), we can find that the system has no equilibrium if ε 1 > 0. The plane ε 1 = 0 excluding the origin in the parameter space is saddle-node bifurcation surface. When ε 1 decrease from this surface, the saddle-node point of Eq (47) becomes a saddle and a node. Then the other bifurcation surfaces are situated in the half space ε 1 < 0. They can be visualized by drawing their trace on the half-sphere when σ > 0 sufficiently small (see Fig 6(a)). We recall that the bifurcation set is a 'cone' based on its trace with S. In Fig 6(b), trace on the S which consists of 3 curves: a curve H om of homoclinic bifurcation, a H of Hopf bifurcation and SN lc of double limit cycle bifurcation. The curve SN lc include two points H 2 and H om2 and the curve SN lc tangent to the curves H and H om . The curves H and H om touch the @S ¼ fðε 1 ; ε 2 ; ε 3 Þjε 1 ¼ 0; ε 2 2 þ ε 2 3 ¼ s 2 g with a first-order tangency at the points BT + and BT − . In the neighbourhood of the BT + and BT − , one can find the unfolding of the cusp-singularity of codimension 2. For system (47), there exists an unstable limit cycle between H and H om near the BT + and an unique stable limit cycle between H and H om near the BT − . In the curved triangle CH 2 H om2 the system has two limit cycles, the inner one unstable and the outer one stable. These two limit cycles coalesce when the ε crosses over the curve SN lc . On the SN lc there exists a unique semistable limit cycle. The more interpretation can be found in literature [24,28].

Bifurcation diagram and simulation
According to analysis and Theorem 0.11, we know that system (5) undergoes Bogdanov-Takens bifurcation of codimension 2. In this section we will choose the parameters β and b as bifurcation parameters to present the bifurcation diagram by simulations. In the proof of Theorem 0.11, we make a series of changes of variables and time, so there will be different situations with different signs of h 0 (I). According to the positive and negative coefficients of XY term in the normal form Eq (47), we denote the Bogdanov-Takens bifurcation of codimension 2 as BT + and BT − respectively.
Taking A = 3, d = 0.3, α = 0.5, μ 0 = 1.5, μ 1 = 3, a = 0, we find that (β, b) = (0.367004, 0.183323) satisfying the conditions in Theorem 0.11, then we use (β, b) to unfold the Bogdanov-Takens bifurcation of codimension 2. By simulation, we obtain the bifurcation diagram in plane (β, b) shown as Fig 7, the blue dash (solid) curve represents saddle-node bifurcation (neutral saddle), the green (blue solid) curve represents supercritical (subcritical) Hopf Dynamic behaviors of a modified SIR model bifurcation and the parameter space (β, b) is divided into differen areas by these bifurcation curves. There are two Bogdanov-Takens bifurcation points, BT − (0.367004, 0.183323) and BT + (0.321298, 0.069049). In order to distinguish these two points, we get some phase diagrams of system when β and b located in different area of (β, b) shown as Fig 8. In Fig 8(a), β, b are located in the area between saddle-node bifurcation and subcritical Hopf bifurcation curve and the epidemic equilibrium E 2 is a unstable focus. In the IV, the phase diagram of system is one of the cases shown as (b), (c) and (h). There is an unstable limit cycle (black curve) near the epidemic equilibrium E 2 in Fig 8(b) and two limit cycles in Fig 8(h) with the inner one unstable and the other one stable. When β and b are located in II, the phase diagram of system is one of the cases as shown in (d) and (e) and there is a stable limit cycle in Fig 8  (e). When β and b are located in I or III, the phase portraits are similar to the cases of (f) and (g), respectively. In the case (f), system (5) has a unique epidemic equilibrium and a stable limit cycle.
In the small neighborhood of BT + , we know that the unstable limit cycle bifurcating from Hopf bifurcation curve disappears from the homoclinic loop, and from   (β, b). There are two types of Bogdanov-Takens bifurcation, BT + and BT − . The green curve represents supercritical Hopf bifurcation, the red curve represents subcritical Hopf bifurcation. The blue dash (solid) curve represents saddle-node bifurcation (neutral saddle curve). https://doi.org/10.1371/journal.pone.0175789.g007 Dynamic behaviors of a modified SIR model  Fig 9, we now know that two limit cycles bifurcating from the semi-stable cycle with one disappearing from the Homoclinic loop and the another disappearing from the Hopf bifurcation curve. We therefore obtain b H 2 > b C > b H om2 by the order these two limit cycle vanishing with different values of b.
Then we obtain the bifurcation diagrams of system (5) near the Bogdanov-Takens bifurcation and the phase portraits in some regions of parameters shown as Figs 10 and 11 respectively.
From the above dynamical analysis, we know that system (5) has complex dynamic behavior even though a = 0. For system (5), we also find the same phenomenon by the simulation as shown in the Fig 12 for the case a 6 ¼ 0. In the simulation, the parameters excluding a are the same as the simulation setting. From Fig 12, we find that the region D 2 and the distance between BT + and BT − becomes small when a becomes lager, which means that choosing a as one other bifurcation parameter can unfold the system (5) and system (5) may undergo the Bogdanov-Takens bifurcation of codimension 3.

Discussion and application
In this paper we consider the SIR model with the nonmonotone incidence rate due to the intervention strategies and nonlinear recovery rate considering the hospitalization conditions. From Theorem 0.1, we know that system (4) undergoes backward bifurcation. In Theorem 0.3, we get the necessary and sufficient condition of backward bifurcation is b < Aðm 1 À m 0 Þ d 2 1 when R 0 ¼ 1, which means that we can eliminate the disease if b > Aðm 1 À m 0 Þ d 2 1 and b < dd 1 A i.e we need enough number of hospital beds. From the Lemma 0.8, we know that if b > m 1 À m 0 À 2d b , system (5) will not have periodic solution, and the endemic equilibrium E 2 is stable. We then discuss Hopf bifurcation and BT bifurcation for system (5) and present in details about these bifurcations in the case a = 0 and present the bifurcation diagrams in Figs 8 and 10.
From the discussion we get Lemma 0.9, which implies that I 2 (b) is a monotone decreasing function of b. Hence, increasing the number of beds can only reduce the number of the total infected individuals, but can not eliminate the diseases as shown in Fig 13(a) if R 0 > 1. If R 0 < 1, from Fig 3 and Eq (36), we know that if b > b c ¼ f À D we can eliminate the disease, and Dynamic behaviors of a modified SIR model  Dynamic behaviors of a modified SIR model In Fig 16, A = 3, β = 0.375, α = 0.5, μ 0 = 0.5 we plot the phase portraits in plane (S, I) with different d, in these cases bA dðdþaþm 0 Þ < 1, and find that there is an unstable limit cycle near the E 2 when d = 1.483783. From the above stimulation, we know that Therorem 0.4 is not ture when À 2 ffiffi ffi a p < c < 0. According to an early SIR model with nonmonotone incidence rate in the literature [19], the dynamics of the system are completely determined by R 0 , which means that the disease will be eliminated if R 0 < 1, otherwise the disease persist. Contrasting to their work and the other results for classic epidemic models, we find that the nonlinear recovery rate is also an important factor which leads to very complicated dynamics. Moreover, we find that R 0 is not enough to determine the dynamic behavior in system (5). By simulations, we predict that there would exist b 1c in system (3)? which has the same role as b c . Hopefully we can explore more relationships between the intervention actions, hospitalization conditions and spread of diseases, to provide the guidelines for public and desicion makers.