Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

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

Dynamic behaviors of a modified SIR model in epidemic diseases using nonlinear incidence and... a1111111111 a1111111111 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 hospi- OPENACCESS talization condition variation in the spread of diseases. By analyzing the existence and sta- bility of the equilibria, we found that the basic reproduction numberR is not a threshold Citation: Li G-H, Zhang Y-X (2017) Dynamic behaviors of a modified SIR model in epidemic parameter, and our model undergoes backward bifurcation when there is limited number of diseases using nonlinear incidence and recovery hospital beds. When the saturated coefficient a is set to zero, it is discovered that the model rates. PLoS ONE 12(4): e0175789. https://doi.org/ undergoes the Saddle-Node bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation 10.1371/journal.pone.0175789 of codimension 2. The bifurcation diagram can further be drawn near the cusp type of the Editor: Gui-Quan Sun, Shanxi University, CHINA Bogdanov-Takens bifurcation of codimension 3 by numerical simulation. We also found a Received: February 14, 2017 critical value of the hospital beds b atR < 1 and sufficiently small a, which suggests that Accepted: March 31, 2017 the disease can be eliminated at the hospitals where the number of beds is larger than b . Published: April 20, 2017 The same dynamic behaviors exist even when a6ˆ 0. Therefore, it can be concluded that a sufficient number of the beds is critical to control the epidemic. Copyright:© 2017 Li, Zhang. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Introduction Data Availability Statement: All relevant data are within the paper and its Supporting Information Since the development of the first dynamic model of smallpox by Bernoulli in 1760, various files. mathematical models have been employed to study infectious diseases [1] in order to reveal Funding: This work is supported by the National the underlying spread mechanisms that influence the transmission and control of these dis- Science Foundation of China (11331009), eases. Among them, Kermack and Mckendrick [2] initiated a famous SIR type of compartmen- Research Project Supported by Fund Program for tal model in 1927 for the plague studies in Mumbai, and succeed in unveiling its the Scientific Activities of Selected Returned epidemiological transition. Since then, mathematical modeling has become an important tool Overseas Professionals in Shanxi Province and the National Sciences Foundation of Shanxi Province to study the transmission and spread of epidemic diseases. (2015011009). 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 mono- Competing interests: The authors have declared that no competing interests exist. tonically 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 Dynamic behaviors of a modified SIR model epidemic models [3]. In those models, the dynamics of models are relatively simple and almost determined by the basic reproduction numberR : the disease will be eliminated ifR < 1, oth- 0 0 erwise, 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 con- trolling 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±15]. Hence, it is essential to expand the modeling studies to the investi- gation 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 −mI a E a I a H 1 2 3 et al [16], Cui et al [17, 18] usedβe , be and c − c f(I) to study the impact of 1 2 media coverage on the dynamics of epidemic models, respectively. However, people may adjust their behaviors according to these government intervention. Therefore, Ruan and Xiao KIS kI S [19] set incidence function in the form of f …I†S ˆ (a special case of [20, 21]) to 2 q 1‡aI 1‡aI 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±24]) started to corpo- rate the limited medical resources into the spread dynamics of infectious diseases. In the litera- ture [22], Wang and Ruan first introduced a piece-wise treatment function in an SIR model, r I > 0; T…I† ˆ 0 I ˆ 0: where the maximal treatment capacity was used to cure infectives so that the epidemic of dis- ease can be controlled. This situation occurs only if the infectious disease needs to be elimi- nated 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 rI 0  I  I ; T…I† ˆ …1† rI I > I : 0 0 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, partic- ularly 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 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 2 / 28 Dynamic behaviors of a modified SIR model infectives. …2† m ˆ m…b; I† ˆ m ‡…m m † 0 1 0 I ‡ b where μ is the minimum per capita recovery rate, and μ the maximum per capita recovery rate. 0 1 They chose the standard incidence rate and discovered the complicated dynamics including Sad- dle-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 bI b…I† ˆ ; aI ‡ cI ‡ 1 p where a is positive constant and c > 2 a (so that aI + 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 hos- pitalization conditions as the literature [24], and the following model is then established, dS bSI > ˆ A dS ; dt 1‡ cI ‡ aI dI bSI …3† ˆ dI aI m…b; I†I; dt 1‡ cI ‡ aI dR ˆ m…b; I†I dR; dt p Fig 1. Graphs of incidence rate function. (a) c 0; (b) 2 a < c < 0. https://doi.org/10.1371/journal.pone.0175789.g001 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 3 / 28 Dynamic behaviors of a modified SIR model where A is the recruitment rate of the susceptible population, d the natural death of the popula- tion, andα the per capita disease-induced death rate, respectively. 3+ 1 For system (3), the cone R is a positive invariant. The C smoothness of the right side of 3+ system (3) implies the local existence and uniqueness of the solution with initial values in R . Adding up the three equations of system (3), we get N (t) = A − dN −αI. Therefore, all solu- tions in the first octant approach, entering or staying inside the set, are defined by D ˆ f…S; I; R†jS  0; I  0; R  0; S‡ I ‡ R  g: 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 c6ˆ 0, but c is in small neighborhood of zero, the behaviors of model still exist. Our model thus becomes, dS bSI ˆ A dS ; dt 1‡ aI dI bSI b ˆ dI aI m ‡…m m † I; …4† 0 1 0 > dt 1‡ aI b‡ I > dR b ˆ m ‡…m m † I dR: 0 1 0 dt b‡ I 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. dS SI ˆ A dS b dt 1‡ aI …5† > dI SI b ˆ b dI aI m ‡…m m † I: 0 1 0 dt 1‡ aI b‡ I We find equilibria by setting the right hand of system (5) equal to zero: SI A dS b ˆ 0; < 2 1‡ aI …6† > SI b dI aI m…b; I†I ˆ 0: 1‡ aI Obviously, a trivial solution of Eq (6) is E …S; I† ˆ ; 0 , a disease free equilibrium(DFE). For E , by using the formula in [25], one can calculate the reproduction number bA R ˆ : …7† d…a‡ d‡ m † For any positive equilibrium E(S, I), also called endemic equilibrium, when exists, its S and I PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 4 / 28 Dynamic behaviors of a modified SIR model coordinates satisfy A …d‡ a‡ m…b; I††I …8† S…I† ˆ ; where the I coordinate will be the positive root of the following cubic equation 3 2 f…I† ˆ A I ‡BI ‡C I ‡D ˆ 0; …9† where A ˆ add ; B ˆ badd ‡ bd ; 0 1 0 C ˆ bbd ‡ dd bA; D ˆ bdd …1 R †; 1 0 1 0 d ˆ d‡ a‡ m ; i ˆ 0; 1: i i Let 2 2 A ˆ B 3AC; B ˆ BC 9AD; C ˆ C 3BD DenoteΔ the discriminant of f(I) with respect to I, then D ˆ B 4A C: 0 1 Note that f…0† ˆ bdd …1 R † and f …0† ˆ C . As shown in Fig 2, we have the following cases 1 0 about the positive roots of f(I): Case 1:R > 1. In this case,D < 0. It is found that there is a unique positive root of f(I) = 0, regardless of the sign of C from Fig 1(c) and 1(d). Case 2:R < 1. In this case,D > 0. IfC > 0, equation f(I) = 0 has no positive solution (see Fig 1(a)). IfC < 0, similar to Lemma 2.1 described by Huang and Ruan [26], the following conclusions can be drawn as shown by Fig 1(b). (a) Δ < 0, there are two positive solutions of the equation. (b) Δ = 0, there is unique positive solution of the equation. (c) Δ > 0, we find that Eq (9) has no positive solution. Case 3:R ˆ 1, Eq (9) becomes 3 2 f…I† ˆ A I ‡BI ‡C I ˆ 0 …10† IfC < 0, Eq (9) has a unique positive root. IfC > 0, Eq (9) has no positive root. Note that d…m m † 1 0 C > 0 means b > . Thus, we get the following theorem about the equilibrium of the bd model. Theorem 0.1. For system (5) with positive parameters, (1) the disease-free equilibrium E always exists, (2) whenR > 1, system has a unique endemic equilibrium, (3) whenR < 1, and (a) C > 0, system (5) does not have endemic equilibrium. PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 5 / 28 Dynamic behaviors of a modified SIR model Fig 2. The positive roots of f(I). (a)R < 1;C > 0; (b)R < 1;C < 0; (c)R > 1;C > 0; (d)R > 1;C < 0. 0 0 0 0 https://doi.org/10.1371/journal.pone.0175789.g002 (b) C < 0,Δ < 0, there exist two endemic equilibria E …I ;S † andE …I ;S †, and 0 1 1 1 2 2 2 whenΔ = 0 these two endemic equilibria coalesce into the same endemic equilibrium E ; otherwise system (5) has no endemic equilibrium. (4) whenR ˆ 1, there exists a unique endemic equilibrium if and only ifC < 0 i.e. d…m m † 1 0 b < ; otherwise there is no endemic equilibrium. bd PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 6 / 28 Dynamic behaviors of a modified SIR model Remark 0.2. By calculation, we get that D …a† ˆ 3b d D‡ O…a†; 0 0 …11† 2 2 2 D ˆ b d b ‡ 2b…2bAd dd d bAd †b‡…dd bA† : 0 0 1 1 0 When a is sufficiently small and tends to zero, the sign ofΔ will be determined by the zero power 2 c of a. Therefore, Δ (a)!Δ (0) = −3β δ Δ as a! 0. In addition, Δ = 0 if and only ifR ˆ R . 0 0 0 0 0 Here R ˆ 1 : 4Bbd…d‡ a‡ m † 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 ) , i, j = 1, 2, then a straightfor- ij 2×2 ward calculation gives bI bS 2abSI j ˆ d ; j ˆ ‡ ; 11 12 2 2 1‡ aI 1‡ aI …1‡ aI † …12† bI bS 2abSI …m m †bI …m m †b 1 0 1 0 j ˆ ; j ˆ d ‡ ‡ : 21 22 0 2 2 2 2 1‡ aI 1‡ aI b‡ I …1‡ aI † …b‡ I† Firstly, we present a theorem about the disease-free equilibrium E (A/d, 0). Theorem 0.3. E is an attracting node ifR < 1, and hyperbolic saddle ifR > 1. When 0 0 d …m m † 1 0 R ˆ 1, E is a saddle-node of codimension 1 if b 6ˆ and attracting semi-hyperbolic node b A d …m m † 1 0 of codimension 2 if b ˆ . b A Proof. For system (5), −d and d …R 1† are two eigenvalues of J(E ). So, E is an attracting 0 0 1 0 node ifR < 1, and unstable ifR > 1. WhenR ˆ 1, the second eigenvalue becomes zero. In 0 0 0 order to analyze the behavior of E , we linearize system (5) and use the transformation of bA X ˆ S‡ I, Y = I, dX ˆ dX ‡ P…X; Y†; dt …13† dY m m b A 0 1 2 ˆ ‡ Y ‡ Q…X; Y†; dt b d where P(X, Y) and Q(X, Y) represent the higher order terms. Obviously, E is a saddle-node if 2 2 d …m m † d …m m † 1 0 1 0 b 6ˆ . Otherwise, i.e., b ˆ , applying the center manifold theorem, system (5) 2 2 b A b A becomes dX > ˆ dX ‡ P…X; Y†; dt …14† dY bAa b Ad > 0 3 ˆ ‡ Y ‡ Q …X; Y†; dt d d where Q (X, Y) represents the higher order term. Thus, E is an attracting semi-hyperbolic 1 0 node of codimension 2. PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 7 / 28 Dynamic behaviors of a modified SIR model Theorem 0.4. If dδ >βA, E is globally asymptotically stable. 0 0 Proof. If dδ >βA, it is obvious thatR < 1 andC > 0. From Theorem 0.1 and 0.3, E is 0 0 the unique attracting node of system (5). In order to prove that the disease free equilibrium E is globally and asymptotically stable, we construct the following Liapunov function: A dS dS V…S; I† ˆ ln ‡ I: …15† d A A It is easy to discover that E ; 0 attains the global minimum of the function V(S, I), so V(S, I) > 0. Along system (5), it turns out: A bAI Vj ˆ 2A dS ‡ …d‡ a‡ m…b; I††I: …16† …5† dS d…1‡ aI † Since μ(b, I) > μ for all I 0, we have 2 3 A …bA dd †I dad I 0 0 Vj  2A dS ‡  0: …17† …5† dS d…1‡ aI † The equality V…S; I† ˆ 0 if and only if at E ; 0 . By Poincare-Bendixson theorem, theorem 0.4 is obvious. Let E…S; I† be any endemic equilibrium, one can verify that its characteristic equation can be written as …18† l tr…J †l‡ det…J † ˆ 0; E E where m m 2aI bI 1 0 tr…J † ˆ d‡ bI …d‡ a‡ m…b; I†† ; 2 2 1‡ aI 1‡ aI …b‡ I† …19† det…J † ˆ …b‡ I†f …I†: E 2 …b‡ I† …1‡ aI † Obviously, the signs of the eigenvalues are determined by f …I† and tr(J ). From Fig 2, we 0 0 know that f …I † < 0; f …I † > 0, so E is a hyperbolic saddle and E is an anti-saddle. E is an 1 2 2 1 2 attracting node or focus, if tr(J ) < 0; E is a weak focus or a center, if tr(J ) = 0; E is a repelling E 2 E 2 node or focus, if tr(J ) > 0. So we obtain the following theorem. Theorem 0.5. For system (5), there are two endemic equilibria E , E whenR < 1 and 1 2 Δ < 0. Then the low endemic equilibrium E is a hyperbolic saddle, and the higher endemic equi- 0 1 librium E is an anti-saddle. WhenR > 1 there is a unique endemic equilibrium, which is an anti-saddle. Bifurcation analysis Backward bifurcation d…m m † 1 0 Theorem 0.6. WhenR ˆ 1, system (5) undergoes backward bifurcation if b < ; and sys- bd d…m m † 1 0 tem (5) undergoes forward bifurcation if b > . bd PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 8 / 28 Dynamic behaviors of a modified SIR model Proof. For convenience of the proof, we suppose that the total number of the population is N(t). System (4) becomes the following system dI b…N I R†I b ˆ dI aI m ‡…m m † I; 0 1 0 dt 1‡ aI b‡ I dR b …20† ˆ m ‡…m m † I dR; 0 1 0 dt b‡ I > dN ˆ A dN aI: dt Let V = (I, R, N) , then the disease-free equilibrium is V ˆ 0; 0; and we can write Eq (20) in vector forms as: …21† V ˆ H…V†…V V †; where 0 1 b…N I R† d a m…b; I† 0 0 B 2 C 1‡ aI B C B C H…V† ˆ : …22† B C m…b; I† d 0 B C @ A a 0 d Then, 0 1 d …R 1† 0 0 1 0 B C B C H…V † ˆ m1 d 0 : …23† B C @ A a 0 d We know that the dominant eigenvalue of H(V ) is zero, ifR ˆ 1. It is well known that we can decompose a neighborhood of the disease-free state into stable manifold W and a center manifold W . 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 is tangential to the 0 c eigenvector V at V . Thus, we can assume that W has the following form: c 0 W ˆ fV ‡ aV ‡ Z…a† : V  Z…a† ˆ 0; a  a  a g; …24† 0 0 0 where V is the dominant left eigenvector of H(V ),α > 0 is a constant, and Z: [−α ,α ]! 0 0 0 0 Ran(H(V )) satisfies: Z…0† ˆ Z…0† ˆ 0: …25† da In other words, W 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 V…t† ˆ V ‡ a…t†V ‡ Z…a…t††; …26† PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 9 / 28 Dynamic behaviors of a modified SIR model since W is an invariant, from Eq (21) we have a _…t†V ‡ Z…a…t†† ˆ V…t† dt ˆ H…V…t††‰V…t† V Š 0 0 ˆ H…V ‡ a…t†V ‡ Z…a…t†††‰a…t†V ‡ Z…a…t††Š; Multiplying both sides of the above equation by V and using the following equations: …27† V  Z…a…t†† ˆ 0; V H…V † ˆ 0; V  V ˆ 1; dt and Z(α) = O(α ) we can get that 0 0 a ˆ V  H…V ‡ aV ‡ Z…a††‰aV ‡ Z…a†Š 0 0 3 ˆ V  H…V ‡ aV †‰aV ‡ Z…a†Š‡ O…a † 0 0 3 ˆ V …H…V ‡ aV † H…V ††‰aV ‡ Z…a†Š‡ O…a †: 0 0 0 0 3 Note that [H(V +αV ) − H(V )] is of orderα, then [H(V +αV ) − H(V )]Z(α) is O(α ) and 0 0 0 0 we get that 0 0 3 a _ ˆ aV ‰H…V ‡ aV † H…V †ŠV ‡ O…a †: …28† 0 0 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: 0 2 3 a _ ˆ V  H V a ‡ O…a †; …29† where dH…V ‡ aV † @H H ˆ j ˆ V j ; …30† aˆ0 VˆV da @V which gives the rate of change of the vector field as the disease invades. Hence, the number h ˆ V H V …31† determines whether the disease can invade whenR ˆ 1, and hence gives the sign of the bifur- cation. For our system, by computation we can get the V and V as follows: m a 1 0 0 1 0 0 V ˆ I ; I ; I ; V ˆ ; 0; 0 ; d d I and 0 0 0 m1 0 a 0 0 0 H ˆ H I ‡ H I H I ; 1 2 3 d d where 0 1 m m 1 0 0 1 0 1 b‡ 0 0 b 0 0 b 0 0 B C B C B C B C 0 B C 0 B C 0 B C m m H ˆ 1 0 ; H ˆ 0 0 0 ; H ˆ 0 0 0 : …32† B C B C B C 1 2 3 0 0 B C @ A @ A @ A 0 0 0 0 0 0 0 0 0 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 10 / 28 Dynamic behaviors of a modified SIR model Then we can get that h ˆ V H V m m b…d‡ a‡ m † 1 0 1 0 ˆ I : …33† b d According to [27], we know that system (5) undergoes backward bifurcation, when h > 0, i.e., d…m m † d…m m † 1 0 1 0 b < ; and system (5) undergoes forward bifurcation, when h < 0, i.e., b > . bd1 bd1 Proposition 0.7. WhenR passes through R and tr(I )6ˆ 0, system (5) with a! 0 under- goes a saddle-node bifurcation. WhenR ˆ R , E is a saddle-node if tr(I )6ˆ 0, and E is a cusp 0 0 if tr(I ) = 0. Proof. As a! 0,Δ (a)!Δ (0).Δ (0) = 0 means thatR ˆ R and the two endemic equilib- 0 0 0 0 0 ria E and E coalesce at E . Two eigenvalues of Jacobian matrix J(E ) are 0 and tr(I ) for system 1 2 (5). 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 b …m m †…bb d† 2 3 0 1 2 < X ˆ X ‡ XO…jYj†‡ O…jYj ;jX; Yj † …b‡ I † jTj …34† Y ˆ tr…I †Y ‡ O…jX; Yj † Where T is the non-singular transformation to diagonalize the linear part. Since b < , E is a saddle-node if tr(I )6ˆ 0. Combined with Theorem 0.1, system (5) undergoes saddle-node bifurcation whenR passes through the critical valueR , 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 numberR ˆ 1 defines a straight line C in the (β, b) plane, dd C : b ˆ : C ˆ 0 also defines one branch of the hyperbolic C (see Fig 3), bA dd C : b ˆ f …b† ˆ : B C bd dd A…m m † 1 1 0 The branch of C intersects with C at the point K ; and withβÐaxis at the B 0 2 A d dd point K ; 0 . It is easily found that f is an increasing convex function ofβ in the first quad- rant. Let dd A…m m † ‡ 1 1 0 C ˆ f…b; b†jb ˆ ; b > g; 0 2 dd A…m m † 1 1 0 C ˆ f…b; b†jb ˆ ; b < g; 0 2 A d ‡ ‡ then C ˆ C [ C [ K . Here C and C are two branches of C joint at point K. 0 0 0 0 0 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 11 / 28 Dynamic behaviors of a modified SIR model Fig 3. The bifurcation curves in (β, b) for system (5) when a6ˆ 0. https://doi.org/10.1371/journal.pone.0175789.g003 Define the curveΔ (β, b) = 0 as C , one can verify that 0 D D …K † ˆ 0; D …K† ˆ 0; 0 0 @b @b j 0 ˆ 0; j ˆ 1: K K @b @b Hence, the curve C is tangent to the curve C at the point K and theβÐaxis at the point K dd dd 0 1 when < b < . A A dd If b ˆ , 2 2 3d …bd ‡ A…d d †† g…b† 0 1 D …b; b† ˆ D …b† ˆ ; 0 0 where 2 2 2 2 2 2 2 g…b† ˆ A a d b 2Aad d b d …4A a…d d † d d †: 1 0 1 0 0 1 0 1 4 3 Denote the discrimination of g(b) = 0 as D ˆ 16A a d d …d d † < 0. Hence the equation 2 0 0 1 A…m m † 1 0 Δ (b) = 0 has a unique real solution, b ˆ , which means that K is the only point at which 0 2 the curve C is tangent to the curve C . If b = 0, D …b; 0† ˆ D …b† ˆ 3d …dd bA† g …b†; 0 0 0 0 1 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 12 / 28 Dynamic behaviors of a modified SIR model where g …b† ˆ d b ‡ 4Aadb 4ad d : 1 0 0 Through computing, we find that equationΔ (β, 0) = 0 has three real solutions q 2 2 2adA 2d A a ‡ ad dd b ˆ ; b ˆ : 0 1;2 A d It is easy to verifyβ >β , so the curve C will not intersect with the abscissa axis when 0 1 dd dd 0 1 b 2 ; . A A For the curve C , note bA dd 12…ad…bA dd †=b‡ bd † …bA dd †…dd bA† 0 0 0 0 1 D b; ˆ bd bd 1 1 …35† 2 2 2 2 2 81a d d …bA dd † …dd bA† 0 0 1 ‡ : 2 2 b d dd dd bA dd 0 1 0 Obviously, if b 2 ; , then D b; > 0. Hence, the curve C is located above the 0 B A A bd dd dd 0 1 curve C for b 2 ; . 0 A A Based on the above the discussion and Theorem 0.1, if we define dd D ˆ f…b; b†jb > ; b > 0g; dd dd A…m m † 0 1 1 0 D ˆ f…b; b†j < b < ; 0 < b < ; D …b; b† < 0g; 2 2 0 A A then there is one endemic equilibria in the region D and two endemic equilibria in the region dd dd D . System (5) undergoes saddle-node bifurcation on the cure C when b 2 ; . The A A backward bifurcation occurs on the C and forward bifurcation occurs on the C . The pitch- 0 0 fork bifurcation occurs when transversally passing through the curve C at the point K. Espe- cially, 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Δ (β, b) = 0, p bA…d d †‡ d …dd bA† 2 bAd …d d †…dd bA† 1 0 0 1 0 1 0 1 …36† b ˆ f …b† ˆ : D 2 bd 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 ) = 0. Note that from Eq (19) we can rewrite tr(J ) as h …I† tr…J † ˆ ; E 2 2 …b‡ I† …1‡ aI † where 4 3 2 2 h …I† ˆ b I ‡ b I ‡ b I ‡ b I ‡ b d; …37† 1 4 3 2 1 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 13 / 28 Dynamic behaviors of a modified SIR model with b ˆ a…3d‡ 2m ‡ 2a†; 4 0 b ˆ a…6bd‡ 3bm ‡ bm ‡ 4ba†‡ b; 3 0 1 2 2 2 b ˆ a…3b d‡ 2b m ‡ 2b a†‡ 2bb‡ d; 2 1 b ˆ b…bb‡ 2d‡ m m †: 1 0 1 m m 2d 1 0 Here, h …I† is a quartic equation. Since b , b and b are non-negative, if b > , 2 3 4 m m 2d 1 0 h (I) > 0, in order to make sure that h (I) = 0 has a positive root, we must have b < . 1 1 This means sufficient number of hospital beds excludes the possibility of the disease oscilla- tion. From the expression of h …I† in Eq (37), it is not an easy task to study the Hopf bifurca- tion 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 3 2 2 h …I† ˆ bI ‡…d‡ 2bb†I ‡ b I ‡ b d: 1 1 One can verify the following lemma m m 2d 1 0 Lemma 0.8. For any positive equilibrium, if b  , we alway have tr(J )<0. In order to study Hopf bifurcation and Bogdanov-Takens bifurcation, we will assume that m m 2d 1 0 b < . Denote the discrimination of h …I† asΔ . Since b < 0, function h …I† must have 1 1 1 1 one negative real root. As shown in Fig 4, It is not difficult to verify function h …I† has two Fig 4. Graph of h(I) with different signs ofΔ when b < 0. When I = H , I = H or I = H = H , Hopf bifurcation occurs. BT bifurcation 1 1 2 m 2 M 2 M m of codimension 2 occurs when I* = H or I* = H and BT bifurcation of codimension 3 occurs when I* = H = H . m M m M https://doi.org/10.1371/journal.pone.0175789.g004 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 14 / 28 Dynamic behaviors of a modified SIR model 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 …I† is determined by the sign of theΔ . When exist, we denote the roots as H and H with H  H . m M M m @I @I 2 2 Lemma 0.9. > 0; 8b > 0 and < 0; 8b > 0. @b @b Proof. From the Eq (9) and the expression of I , direct calculation leads to @I @f @f dd I ‡ bdd 2 0 1 p ˆ = ˆ > 0; @b @b @I b D @I 1 @C 1 @C @D ˆ ‡p C 2B : @b 2B @b @b @b @I @C One can find that ˆ bd > 0. We will prove < 0 in two cases. IfR < 1, then @b 1 @b 0 @D ˆ dd bA > 0. Recall the analysis of the existence of the equilibria we know that the @b @I @I 2 2 C < 0, so the < 0. IfR > 1, then lim ˆ 0 and 0 b!‡1 @b @b " # 3=2 2 @I D @D @ D 2 0 0 0 ˆ 2D 2 0 2 8B @b @ b @ b ˆ 2b A…m m †…bA dd † > 0; 1 0 1 @I …b† so the is an increasing function of b with supremum 0 in the (0, +1), so for8b > 0 @b @I …b† < 0. @b Theorem 0.10. For system (5) with a = 0, generic Hopf bifurcation could occur if I = H , 2 m I = H or I = H = H . 2 M 2 m M Proof. We only need to verify the transversality condition. Letγ = tr(I )/2 be the real part of the two solutions of Eq (18), when a = 0. If I = H or I = H = H , we considerβ as the bifurcation parameter and fix all other 2 M 2 M m parameters. Then dg 1 @tr…I …b†; b† @I …b† @tr…I …b†; b† 2 2 2 ˆ ‡ ; db 2 @I @b @b ^ ^ bˆb 2 bˆb @tr…I …b†; b† 2h…H † 2 M ˆ < 0; @I ^ …b‡ H † 2 bˆb 3 2 2 @tr…I …b†; b† H ‡ 2bH ‡ b H 2 M M M ˆ < 0: @b ^ …b‡ H † bˆb @I dg From Lemma 0.9, we have j > 0, so < 0. bˆb @b db If I = H or I = H = H , we consider b as the bifurcation parameter and fix all the other 2 m 2 M m parameters. Then dg 1 @tr…I …b†; b†@I …b† @tr…I …b†; b† 2 2 2 ˆ ‡ ; db 2 @I @b @b ^ ^ bˆb 2 bˆb @tr…I …b†; b† 2h…H † 2 m ˆ > 0; @b ^ …b‡ H † bˆb m @tr…I …b†; b† …b H †…m m † 2 m 0 1 ˆ < 0: @b ^ …b‡ H † bˆb m PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 15 / 28 Dynamic behaviors of a modified SIR model Fig 5. Graphs of Bifurcation curve in parameters plane (β, b) and the phase trajectory for system (5). (a) Curve q(β, b) = 0. The green curve is supercritical Hopf bifurcation; The red curve is subcritical Hopf bifurcation.σ becomes 0 at the DH point. (b) Two limit cycles bifurcation from the weak focus E . https://doi.org/10.1371/journal.pone.0175789.g005 @I dg From Lemma 0.9, we have j < 0, so < 0 (one can verify that H < b). Then the proof @b bˆb dm 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 . If E is a weak focus, Hopf bifurcation can happen, otherwise system does not undergo 2 2 Hopf bifurcation. Because system (5) is analytic when a = 0, E can only be weak focus or cen- ter. We can distinguish these two types of singularities by calculating Lyapunov coefficients and verifying it through numerical simulation. Taking the resultant of f(I) and h …I† with respect to I, we can get the curve q(β, b) in the space (β, b), and plot the algebraic curve q(β, b) = 0 by fixing other parameters A, d, μ ,α and μ . Choose A = 3, d = 0.3,α = 0.5, μ = 1.5, μ = 3 and plot q(β, b) = 0 in the plane (β, b) as 0 0 1 shown in Fig 5(a). The green curve (δ < 0) represents supercritical Hopf bifurcation; the red curve corresponding toδ > 0 represents subcritical Hopf bifurcation. 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 (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 . Bogdanov-Takens bifurcation From Theorem 0.1 we know that the two equilibria E and E coalesce at the equilibria 1 2 E (S , I ) whenR ˆ R , if a = 0, where 0 0 S ˆ ; d‡ bI d…d‡ a‡ m † bA‡ bb…d‡ a‡ m † 0 1 I ˆ : 2b…d‡ a‡ m † We can find that det(I ) = 0 in Eq (18) ifR ˆ R . From Proposition 0.7, we know that E is a 0 0 saddle-node point if tr(I )6ˆ 0. If tr(I ) = 0, Eq (18) has a zero eigenvalue with multiplicity 2, PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 16 / 28 Dynamic behaviors of a modified SIR model 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 (I )6ˆ 0, then E is a Bogdanov-Takens point of codimension 2, and system (5) localized at E is topologically equivalent to X ˆ Y; …38† : 3 Y ˆ X ‡ Sign…h…I ††XY ‡ O…jX; Yj †: Proof. Changing the variables as x = S − S , y = I − I , then system (5) becomes dx > ˆ …d‡ bI †x bS y bxy; dt …39† dy …m m †bI 1 0 2 : ˆ bI x‡ y‡ Cy ‡ bxy‡ O…jx; yj †; dt …b‡ I † where …m m †b …m m †bI 1 0 1 0 C ˆ : …40† 2 3 …b‡ I † …b‡ I † tr(I ) = 0 and det(I ) = 0 imply that …m m †bI 2 1 0 b S I …d‡ bI † ˆ 0; d‡ bI ˆ ; …41† …b‡ I † so the generalized eigenvectors corresponding toλ = 0 of Jacobian matrix in system (5) at the point E are 0 0 …42† V ˆ … d bI ; bI † ; V ˆ …1; 0† : 1 2 Let T = (T ) = (V , V ), then under the non-singular linear transformation ij 2×2 1 2 ! ! x X ˆ T ; y Y where |T| = −βI < 0. System (39) becomes X ˆ Y ‡ a X ‡ bXY; …43† : 3 Y ˆ a X ‡ dbXY ‡ O…jX; Yj †; here …m m †b …m m †bI 2 1 0 1 0 2 …d‡ bI †b I b I 2 3 …b‡ I † …b‡ I † a ˆ ; bI bI …bI ‡ d†…bb d† a ˆ : b‡ I PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 17 / 28 Dynamic behaviors of a modified SIR model Using the near-identity transformation a X 12 2 …44† u ˆ X ; v ˆ Y ‡ a X ; and rewrite u, v into X, Y, we obtain X ˆ Y ‡ O…jX; Yj †; …45† : 3 Y ˆ M X ‡ M XY ‡ O…jX; Yj †: 21 22 To consider the sign of M , note that bI …bI ‡ d†…bb d† M ˆ a ˆ : 21 21 b‡ I d…m m † 1 0 For system (5), the condition of the existence of endemic equilibrium is b < , hence, bd M < 0. Then we will determine the sign of M by 21 22 M ˆ a ‡ 2a 22 22 11 2 …m m †b …m m †bI 2 1 0 1 0 …d‡ bI †b I b I 2 3 …b‡I † …b‡I † ˆ db 2 bI bh…I † ˆ : …b‡ I † If h (I )6ˆ 0, we make a change of coordinates and time and preserve the orientation by time M M a 21 21 22 X ! X; Y ! Y; t ! j jt …46† a a M 22 21 then system (5) is topologically equivalent to the normal form Eq (38). From Theorem 0.11, we know that if a = 0, endemic equilibrium E is a Bogdanov-Takens 0  0 point of codimension 2 whenΔ = 0, h(I ) = 0 and h (I )6ˆ 0. If h (I ) = 0, E may be a cusp of codimension 3. In [28], a generic unfolding with the parametersε = (ε ,ε ,ε ) of the codimension 3 cusp 1 2 3 singularity is C equivalent to X ˆ Y; …47† : 4 2 3 Y ˆ ε ‡ε Y ‡ε XY ‡ X X Y ‡ O…jX; Yj †: 1 2 3 About system (47), we can find that the system has no equilibrium ifε > 0. The planeε = 1 1 0 excluding the origin in the parameter space is saddle-node bifurcation surface. Whenε 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ε < 0. They can be visualized by drawing their trace on the half-sphere 2 2 2 2 S ˆ f…ε ;ε ;ε †jε < 0;ε ‡ε ‡ε ˆ s g; …48† 1 2 3 1 1 2 3 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 of homoclinic bifurcation, om a H of Hopf bifurcation and SN of double limit cycle bifurcation. The curve SN include two lc lc PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 18 / 28 Dynamic behaviors of a modified SIR model Fig 6. The bifurcation diagram of BT of codimension 3. (a) The parameter space and the trace of the bifurcation diagram on the S(  0); (b) The sign of the BT is positive if the coefficient of the term XY in the norm form is positive, otherwise it is negative [24]. https://doi.org/10.1371/journal.pone.0175789.g006 points H and H and the curve SN tangent to the curves H and H . The curves H and 2 om2 lc om 2 2 2 H touch the @S ˆ f…ε ;ε ;ε †jε ˆ 0;ε ‡ε ˆ s g with a first-order tangency at the om 1 2 3 1 2 3 + − + − 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 near the BT and an unique stable limit cycle between H and H near the om om BT . In the curved triangle CH H the system has two limit cycles, the inner one unstable 2 om2 and the outer one stable. These two limit cycles coalesce when theε crosses over the curve SN . On the SN there exists a unique semistable limit cycle. The more interpretation can be lc lc 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 The- orem 0.11, we make a series of changes of variables and time, so there will be different situa- tions with different signs of h (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, μ = 1.5, μ = 3, a = 0, we find that (β, b) = (0.367004, 0 1 0.183323) satisfying the conditions in Theorem 0.11, then we use (β, b) to unfold the Bogda- nov-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 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 19 / 28 Dynamic behaviors of a modified SIR model + − Fig 7. The bifurcation diagram in plane (β, 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 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 is a unstable focus. In the IV, the phase dia- gram 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 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 dia- gram 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 Fig 8(b) and 8(c), we can observe that the homoclinic loops are located in IV. Otherwise, from Fig 8(d) and 8(e), we can obtain that the homoclinic loops are located in II which is in the small neighborhood of BT . Hence, the Hopf bifurcation curve and homoclinic loops switch their positions at some point C. In order to figure out the relative positions of C, H and Hom , as shown in Fig 9 we 2 2 change the value ofβ and plot the bifurcation diagram on (β, b) with different b. PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 20 / 28 Dynamic behaviors of a modified SIR model Fig 8. The phase diagram of system (5). The blue curve represents unstable manifold, green curve represents stable manifold. (a) b = 0.1,β = 0.339; (b) b = 0.1,β = 0.340; (c) b = 0.1,β = 0.3415; (d) b = 0.18,β = 0.367; (e) b = 0.18,β = 0.3737; (f) b = 0.21,β = 0.3815; (g) b = 0.21,β = 0.4; (h) b = 0.1587,β = 0.3683. In the (b), there is an unstable limit cycle marked black curve near the epidemic equilibrium E . In the (e) and (f), there is a stable limit cycle marked red curve. In the (h), we find that there are two limit cycle, the small one is unstable, the another one is stable. https://doi.org/10.1371/journal.pone.0175789.g008 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 21 / 28 Dynamic behaviors of a modified SIR model Fig 9. Bifurcation diagram in (β, I) with different b. The red dash(solid) represents unstable epidemic equilibrium(limit cycle). The blue curve represents stable epidemic equilibrium or limit cycle. (a) b = 0.145; (b) b = 0.14. https://doi.org/10.1371/journal.pone.0175789.g009 By simulations, we get the case (a) and (b) in Fig 9, we now know that two limit cycles bifur- cating 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 > b > b H C H 2 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 bifurca- tion 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 behav- ior 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 a6ˆ 0. In the simulation, the parameters excluding a are the same as the simulation setting. From Fig 12, we find that the region D 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 A…m m † 1 0 0.3, we get the necessary and sufficient condition of backward bifurcation is b < when A…m m † dd 1 0 1 R ˆ 1, which means that we can eliminate the disease if b > and b < i.e we need m m 2d 1 0 enough number of hospital beds. From the Lemma 0.8, we know that if b > , system (5) will not have periodic solution, and the endemic equilibrium E is stable. We then discuss Hopf bifurcation and BT bifurcation for system (5) and present in details about these bifurca- tions 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 (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) ifR > 1. If R < 1, from Fig 3 and Eq (36), we know that if b > b ˆ f we can eliminate the disease, and 0 c D PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 22 / 28 Dynamic behaviors of a modified SIR model Fig 10. Bifurcation digram near the Bogdanov-Takens. https://doi.org/10.1371/journal.pone.0175789.g010 Fig 11. Phase portraits for parameters in different regions of Fig 10. https://doi.org/10.1371/journal.pone.0175789.g011 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 23 / 28 Dynamic behaviors of a modified SIR model Fig 12. The curve q(β, b) = 0 with different values of a. The blue curve, yellow curve, red curve and green curve are drawn according to a = 0, a = 0.5, a = 1 and a = 2 respectively. https://doi.org/10.1371/journal.pone.0175789.g012 these rich dynamics finally disappear through the saddle-node bifurcation when b = b as shown in Figs 7 and 13(b) and 13(c). For system (3) with a6ˆ 0, taking A = 3, d = 0.3,β = 0.5, a = 0.2, μ = 3.1728627 and μ = 1.5 1 0 we get the bifurcation diagram with different values for c as shown as Figs 14 and 15. In Fig 14, the types of BT-bifurcation are the same, however, there are two types of BT bifurcations in the Fig 15. Fig 13. Bifurcation in the plane (b, I) with differentμ .β = 0.39, a = 0. (a)R ˆ 1:02; (b)R ˆ 0:98; (c)R ˆ 0:95. 1 0 0 0 https://doi.org/10.1371/journal.pone.0175789.g013 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 24 / 28 Dynamic behaviors of a modified SIR model Fig 14. The bifurcation diagram of system (3) in parameters plane (β, b). A = 3, d = 0.3,β = 0.5, c = 0.185, a = 0.2,μ = 3.1728627,μ = 1 0 ‡ ‡ 1.5,BT …0:353073; 0:0925301†;BT …0:375; 0:137406†. 1 2 https://doi.org/10.1371/journal.pone.0175789.g014 Fig 15. The bifurcation diagram of system (3) in parameters plane (β, b). A = 3, d = 0.3,β = 0.5, c = 0.1, a = 0.2,μ = 3.1728627,μ = 1 0 + − 1.5, BT (0.337066, 0.072821), BT (0.382572, 0.172627). https://doi.org/10.1371/journal.pone.0175789.g015 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 25 / 28 Dynamic behaviors of a modified SIR model Fig 16. Phase portraits of system (3) in the plane (S, I) with different d. A = 3,β = 0.375,α = 0.5,μ = 0.5,μ = 0.7, a = 0.7, b = 0.01, c = 0 1 −1.65. (a) d = 1.49; (b) d = 1.483783; (c) d = 1.4. In case (a) (b) and (c), E is always a saddle and E is a stable node. E is a stable node in 1 0 2 case (b) and (c), while it is a unstable node in case (a). There is an unstable limit cycle near E in case (b). https://doi.org/10.1371/journal.pone.0175789.g016 In Fig 16, A = 3,β = 0.375,α = 0.5, μ = 0.5 we plot the phase portraits in plane (S, I) with bA different d, in these cases < 1, and find that there is an unstable limit cycle near the E d…d‡a‡m † when d = 1.483783. From the above stimulation, we know that Therorem 0.4 is not ture when p 2 a < c < 0. According to an early SIR model with nonmonotone incidence rate in the literature [19], the dynamics of the system are completely determined byR , which means that the disease will be eliminated ifR < 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 thatR is not enough to determine the dynamic behavior in system (5). By simulations, we predict that there would exist b in system (3)? which has the same role as b . Hopefully we can explore more 1c c relationships between the intervention actions, hospitalization conditions and spread of dis- eases, to provide the guidelines for public and desicion makers. Acknowledgments The authors would like to thank anonymous referees for very helpful suggestions and com- ments which led to improvements of our original manuscript. Author Contributions Conceptualization: GHL YXZ. Formal analysis: GHL YXZ. Writing ± original draft: GHL YXZ. Writing ± review & editing: GHL YXZ. References 1. Brauer F, Castillo-Chavez C. Mathematical models in population biology and epidemiology[M]. New York: Springer, 2001. 2. Kermack W O, McKendrick A G. A contribution to the mathematical theory of epidemics[C] Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Soci- ety, 1927, 115(772): 700±721. PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 26 / 28 Dynamic behaviors of a modified SIR model 3. Hethcote H W. The mathematics of infectious diseases[J]. SIAM review, 2000, 42(4): 599±653. https:// doi.org/10.1137/S0036144500371907 4. Gumel A B, Ruan S, Day T, et al. Modelling strategies for controlling SARS outbreaks[J]. Proceedings of the Royal Society of London B: Biological Sciences, 2004, 271(1554): 2223±2232. https://doi.org/10. 1098/rspb.2004.2800 5. Li Li. Monthly periodic outbreak of hemorrhagic fever with renal syndrome in China[J]. Journal of Biologi- cal Systems, 2016, 24, 519±533. https://doi.org/10.1142/S0218339016500261 6. Li Li. Patch invasion in a spatial epidemic model[J]. Applied Mathematics and Computation, 2015, 258, 342±349. https://doi.org/10.1016/j.amc.2015.02.006 7. Gui-Quan Sun, et al, Transmission dynamics of cholera: Mathematical modeling and control strategies [J]. Commun. Nonlinear Sci. Numer. Simulat., 2017, 45, 235±244. https://doi.org/10.1016/j.cnsns. 2016.10.007 8. Ming-Tao Li, Jin Z., Sun G., Zhang J.. Modeling direct and indirect disease transmission using multi- group model[J]. J. Math. Anal. Appl., 2017, 446, 1292±1309. https://doi.org/10.1016/j.jmaa.2016.09. 9. Ming-Tao Li, Sun G., Wu Y., Zhang J. and Jin Z.. Transmission dynamics of a multi-group brucellosis model with mixed cross infection in public farm[J]. Applied Mathematics and Computation, 2014, 237, 582±594. https://doi.org/10.1016/j.amc.2014.03.094 10. Gui-Quan Sun, Zhang Z., Global stability for a sheep brucellosis model with immigration[J], Applied Mathematics and Computation, 2014, 246, 336±345. https://doi.org/10.1016/j.amc.2014.08.028 11. Zhi-Qiang Xia, et al, Modeling the transmission dynamics of Ebola virus disease in Liberia[J]. Scientific Reports, 2015, 5, 13857. https://doi.org/10.1038/srep13857 12. Yan-Fang Wu, Li M., Sun G.. Asymptotic analysis of schistosomiasis persistence in models with general functions[J]. Journal of the Franklin Institute, 2016, 353, 4772±4784. https://doi.org/10.1016/j.jfranklin. 2016.09.012 13. Zhang J, Lou J, Ma Z, et al. A compartmental model for the analysis of SARS transmission patterns and outbreak control measures in China[J]. Applied Mathematics and Computation, 2005, 162(2): 909± 924. https://doi.org/10.1016/j.amc.2003.12.131 14. Tang S, Xiao Y, Yuan L, et al. Campus quarantine (Fengxiao) for curbing emergent infectious diseases: lessons from mitigating A/H1N1 in Xi'an, China[J]. Journal of theoretical biology, 2012, 295: 47±58. https://doi.org/10.1016/j.jtbi.2011.10.035 PMID: 22079943 15. Xiao Y, Tang S, Wu J. Media impact switching surface during an infectious disease outbreak[J]. Scien- tific reports, 2015, 5. 16. Liu R, Wu J, Zhu H. Media/psychological impact on multiple outbreaks of emerging infectious diseases [J]. Computational and Mathematical Methods in Medicine, 2007, 8(3): 153±164. https://doi.org/10. 1080/17486700701425870 17. Cui J, Sun Y, Zhu H. The impact of media on the control of infectious diseases[J]. Journal of Dynamics and Differential Equations, 2008, 20(1): 31±53. https://doi.org/10.1007/s10884-007-9075-0 18. Cui J, Tao X, Zhu H. An SIS infection model incorporating media coverage[J]. Rocky Mountain J. Math, 2008, 38(5): 1323±1334. https://doi.org/10.1216/RMJ-2008-38-5-1323 19. Xiao D, Ruan S. Global analysis of an epidemic model with nonmonotone incidence rate[J]. Mathemati- cal biosciences, 2007, 208(2): 419±429. https://doi.org/10.1016/j.mbs.2006.09.025 PMID: 17303186 20. Liu W, Levin S A, Iwasa Y. Influence of nonlinear incidence rates upon the behavior of SIRS epidemio- logical models[J]. Journal of mathematical biology, 1986, 23(2): 187±204. https://doi.org/10.1007/ BF00276956 PMID: 3958634 21. Liu W, Hethcote H W, Levin S A. Dynamical behavior of epidemiological models with nonlinear inci- dence rates[J]. Journal of mathematical biology, 1987, 25(4): 359±380. https://doi.org/10.1007/ BF00277162 PMID: 3668394 22. Wang W, Ruan S. Bifurcations in an epidemic model with constant removal rate of the infectives[J]. Journal of Mathematical Analysis and Applications, 2004, 291(2): 775±793. https://doi.org/10.1016/j. jmaa.2003.11.043 23. Wang W. Backward bifurcation of an epidemic model with treatment[J]. Mathematical biosciences, 2006, 201(1): 58±71. https://doi.org/10.1016/j.mbs.2005.12.022 PMID: 16466756 24. Shan C, Zhu H. Bifurcations and complex dynamics of an SIR model with the impact of the number of hospital beds[J]. Journal of Differential Equations, 2014, 257(5): 1662±1688. https://doi.org/10.1016/j. jde.2014.05.030 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 27 / 28 Dynamic behaviors of a modified SIR model 25. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[J]. Mathematical biosciences, 2002, 180(1): 29±48. https://doi.org/10.1016/S0025-5564(02)00108-6 PMID: 12387915 26. Huang J, Ruan S, Song J. Bifurcations in a predator?Cprey system of Leslie type with generalized Hol- ling type III functional response[J]. Journal of Differential Equations, 2014, 257(6): 1721±1752. https:// doi.org/10.1016/j.jde.2014.04.024 27. Dushoff F, Huang W, Castillo-Chavez C. Backward bifurcation and catastrophe in simple models of fatal disease[J]. Journal of mathematical biology, 1998, 36(3): 227±248. https://doi.org/10.1007/ s002850050099 PMID: 9528119 28. Dumortier F, Roussarie R, Sotomayor J. Generic 3-parameter families of vector fields on the plane, unfolding a singularity with nilpotent linear part. The cusp case of codimension 3[J]. Ergodic theory and dynamical systems, 1987, 7(03): 375±413. https://doi.org/10.1017/S0143385700004119 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 28 / 28 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png PLoS ONE Pubmed Central

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

PLoS ONE , Volume 12 (4) – Apr 20, 2017

Loading next page...
 
/lp/pubmed-central/dynamic-behaviors-of-a-modified-sir-model-in-epidemic-diseases-using-0Uyj8Yco2K

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Pubmed Central
Copyright
© 2017 Li, Zhang
ISSN
1932-6203
eISSN
1932-6203
DOI
10.1371/journal.pone.0175789
Publisher site
See Article on Publisher Site

Abstract

a1111111111 a1111111111 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 hospi- OPENACCESS talization condition variation in the spread of diseases. By analyzing the existence and sta- bility of the equilibria, we found that the basic reproduction numberR is not a threshold Citation: Li G-H, Zhang Y-X (2017) Dynamic behaviors of a modified SIR model in epidemic parameter, and our model undergoes backward bifurcation when there is limited number of diseases using nonlinear incidence and recovery hospital beds. When the saturated coefficient a is set to zero, it is discovered that the model rates. PLoS ONE 12(4): e0175789. https://doi.org/ undergoes the Saddle-Node bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation 10.1371/journal.pone.0175789 of codimension 2. The bifurcation diagram can further be drawn near the cusp type of the Editor: Gui-Quan Sun, Shanxi University, CHINA Bogdanov-Takens bifurcation of codimension 3 by numerical simulation. We also found a Received: February 14, 2017 critical value of the hospital beds b atR < 1 and sufficiently small a, which suggests that Accepted: March 31, 2017 the disease can be eliminated at the hospitals where the number of beds is larger than b . Published: April 20, 2017 The same dynamic behaviors exist even when a6ˆ 0. Therefore, it can be concluded that a sufficient number of the beds is critical to control the epidemic. Copyright:© 2017 Li, Zhang. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Introduction Data Availability Statement: All relevant data are within the paper and its Supporting Information Since the development of the first dynamic model of smallpox by Bernoulli in 1760, various files. mathematical models have been employed to study infectious diseases [1] in order to reveal Funding: This work is supported by the National the underlying spread mechanisms that influence the transmission and control of these dis- Science Foundation of China (11331009), eases. Among them, Kermack and Mckendrick [2] initiated a famous SIR type of compartmen- Research Project Supported by Fund Program for tal model in 1927 for the plague studies in Mumbai, and succeed in unveiling its the Scientific Activities of Selected Returned epidemiological transition. Since then, mathematical modeling has become an important tool Overseas Professionals in Shanxi Province and the National Sciences Foundation of Shanxi Province to study the transmission and spread of epidemic diseases. (2015011009). 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 mono- Competing interests: The authors have declared that no competing interests exist. tonically 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 Dynamic behaviors of a modified SIR model epidemic models [3]. In those models, the dynamics of models are relatively simple and almost determined by the basic reproduction numberR : the disease will be eliminated ifR < 1, oth- 0 0 erwise, 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 con- trolling 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±15]. Hence, it is essential to expand the modeling studies to the investi- gation 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 −mI a E a I a H 1 2 3 et al [16], Cui et al [17, 18] usedβe , be and c − c f(I) to study the impact of 1 2 media coverage on the dynamics of epidemic models, respectively. However, people may adjust their behaviors according to these government intervention. Therefore, Ruan and Xiao KIS kI S [19] set incidence function in the form of f …I†S ˆ (a special case of [20, 21]) to 2 q 1‡aI 1‡aI 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±24]) started to corpo- rate the limited medical resources into the spread dynamics of infectious diseases. In the litera- ture [22], Wang and Ruan first introduced a piece-wise treatment function in an SIR model, r I > 0; T…I† ˆ 0 I ˆ 0: where the maximal treatment capacity was used to cure infectives so that the epidemic of dis- ease can be controlled. This situation occurs only if the infectious disease needs to be elimi- nated 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 rI 0  I  I ; T…I† ˆ …1† rI I > I : 0 0 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, partic- ularly 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 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 2 / 28 Dynamic behaviors of a modified SIR model infectives. …2† m ˆ m…b; I† ˆ m ‡…m m † 0 1 0 I ‡ b where μ is the minimum per capita recovery rate, and μ the maximum per capita recovery rate. 0 1 They chose the standard incidence rate and discovered the complicated dynamics including Sad- dle-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 bI b…I† ˆ ; aI ‡ cI ‡ 1 p where a is positive constant and c > 2 a (so that aI + 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 hos- pitalization conditions as the literature [24], and the following model is then established, dS bSI > ˆ A dS ; dt 1‡ cI ‡ aI dI bSI …3† ˆ dI aI m…b; I†I; dt 1‡ cI ‡ aI dR ˆ m…b; I†I dR; dt p Fig 1. Graphs of incidence rate function. (a) c 0; (b) 2 a < c < 0. https://doi.org/10.1371/journal.pone.0175789.g001 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 3 / 28 Dynamic behaviors of a modified SIR model where A is the recruitment rate of the susceptible population, d the natural death of the popula- tion, andα the per capita disease-induced death rate, respectively. 3+ 1 For system (3), the cone R is a positive invariant. The C smoothness of the right side of 3+ system (3) implies the local existence and uniqueness of the solution with initial values in R . Adding up the three equations of system (3), we get N (t) = A − dN −αI. Therefore, all solu- tions in the first octant approach, entering or staying inside the set, are defined by D ˆ f…S; I; R†jS  0; I  0; R  0; S‡ I ‡ R  g: 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 c6ˆ 0, but c is in small neighborhood of zero, the behaviors of model still exist. Our model thus becomes, dS bSI ˆ A dS ; dt 1‡ aI dI bSI b ˆ dI aI m ‡…m m † I; …4† 0 1 0 > dt 1‡ aI b‡ I > dR b ˆ m ‡…m m † I dR: 0 1 0 dt b‡ I 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. dS SI ˆ A dS b dt 1‡ aI …5† > dI SI b ˆ b dI aI m ‡…m m † I: 0 1 0 dt 1‡ aI b‡ I We find equilibria by setting the right hand of system (5) equal to zero: SI A dS b ˆ 0; < 2 1‡ aI …6† > SI b dI aI m…b; I†I ˆ 0: 1‡ aI Obviously, a trivial solution of Eq (6) is E …S; I† ˆ ; 0 , a disease free equilibrium(DFE). For E , by using the formula in [25], one can calculate the reproduction number bA R ˆ : …7† d…a‡ d‡ m † For any positive equilibrium E(S, I), also called endemic equilibrium, when exists, its S and I PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 4 / 28 Dynamic behaviors of a modified SIR model coordinates satisfy A …d‡ a‡ m…b; I††I …8† S…I† ˆ ; where the I coordinate will be the positive root of the following cubic equation 3 2 f…I† ˆ A I ‡BI ‡C I ‡D ˆ 0; …9† where A ˆ add ; B ˆ badd ‡ bd ; 0 1 0 C ˆ bbd ‡ dd bA; D ˆ bdd …1 R †; 1 0 1 0 d ˆ d‡ a‡ m ; i ˆ 0; 1: i i Let 2 2 A ˆ B 3AC; B ˆ BC 9AD; C ˆ C 3BD DenoteΔ the discriminant of f(I) with respect to I, then D ˆ B 4A C: 0 1 Note that f…0† ˆ bdd …1 R † and f …0† ˆ C . As shown in Fig 2, we have the following cases 1 0 about the positive roots of f(I): Case 1:R > 1. In this case,D < 0. It is found that there is a unique positive root of f(I) = 0, regardless of the sign of C from Fig 1(c) and 1(d). Case 2:R < 1. In this case,D > 0. IfC > 0, equation f(I) = 0 has no positive solution (see Fig 1(a)). IfC < 0, similar to Lemma 2.1 described by Huang and Ruan [26], the following conclusions can be drawn as shown by Fig 1(b). (a) Δ < 0, there are two positive solutions of the equation. (b) Δ = 0, there is unique positive solution of the equation. (c) Δ > 0, we find that Eq (9) has no positive solution. Case 3:R ˆ 1, Eq (9) becomes 3 2 f…I† ˆ A I ‡BI ‡C I ˆ 0 …10† IfC < 0, Eq (9) has a unique positive root. IfC > 0, Eq (9) has no positive root. Note that d…m m † 1 0 C > 0 means b > . Thus, we get the following theorem about the equilibrium of the bd model. Theorem 0.1. For system (5) with positive parameters, (1) the disease-free equilibrium E always exists, (2) whenR > 1, system has a unique endemic equilibrium, (3) whenR < 1, and (a) C > 0, system (5) does not have endemic equilibrium. PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 5 / 28 Dynamic behaviors of a modified SIR model Fig 2. The positive roots of f(I). (a)R < 1;C > 0; (b)R < 1;C < 0; (c)R > 1;C > 0; (d)R > 1;C < 0. 0 0 0 0 https://doi.org/10.1371/journal.pone.0175789.g002 (b) C < 0,Δ < 0, there exist two endemic equilibria E …I ;S † andE …I ;S †, and 0 1 1 1 2 2 2 whenΔ = 0 these two endemic equilibria coalesce into the same endemic equilibrium E ; otherwise system (5) has no endemic equilibrium. (4) whenR ˆ 1, there exists a unique endemic equilibrium if and only ifC < 0 i.e. d…m m † 1 0 b < ; otherwise there is no endemic equilibrium. bd PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 6 / 28 Dynamic behaviors of a modified SIR model Remark 0.2. By calculation, we get that D …a† ˆ 3b d D‡ O…a†; 0 0 …11† 2 2 2 D ˆ b d b ‡ 2b…2bAd dd d bAd †b‡…dd bA† : 0 0 1 1 0 When a is sufficiently small and tends to zero, the sign ofΔ will be determined by the zero power 2 c of a. Therefore, Δ (a)!Δ (0) = −3β δ Δ as a! 0. In addition, Δ = 0 if and only ifR ˆ R . 0 0 0 0 0 Here R ˆ 1 : 4Bbd…d‡ a‡ m † 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 ) , i, j = 1, 2, then a straightfor- ij 2×2 ward calculation gives bI bS 2abSI j ˆ d ; j ˆ ‡ ; 11 12 2 2 1‡ aI 1‡ aI …1‡ aI † …12† bI bS 2abSI …m m †bI …m m †b 1 0 1 0 j ˆ ; j ˆ d ‡ ‡ : 21 22 0 2 2 2 2 1‡ aI 1‡ aI b‡ I …1‡ aI † …b‡ I† Firstly, we present a theorem about the disease-free equilibrium E (A/d, 0). Theorem 0.3. E is an attracting node ifR < 1, and hyperbolic saddle ifR > 1. When 0 0 d …m m † 1 0 R ˆ 1, E is a saddle-node of codimension 1 if b 6ˆ and attracting semi-hyperbolic node b A d …m m † 1 0 of codimension 2 if b ˆ . b A Proof. For system (5), −d and d …R 1† are two eigenvalues of J(E ). So, E is an attracting 0 0 1 0 node ifR < 1, and unstable ifR > 1. WhenR ˆ 1, the second eigenvalue becomes zero. In 0 0 0 order to analyze the behavior of E , we linearize system (5) and use the transformation of bA X ˆ S‡ I, Y = I, dX ˆ dX ‡ P…X; Y†; dt …13† dY m m b A 0 1 2 ˆ ‡ Y ‡ Q…X; Y†; dt b d where P(X, Y) and Q(X, Y) represent the higher order terms. Obviously, E is a saddle-node if 2 2 d …m m † d …m m † 1 0 1 0 b 6ˆ . Otherwise, i.e., b ˆ , applying the center manifold theorem, system (5) 2 2 b A b A becomes dX > ˆ dX ‡ P…X; Y†; dt …14† dY bAa b Ad > 0 3 ˆ ‡ Y ‡ Q …X; Y†; dt d d where Q (X, Y) represents the higher order term. Thus, E is an attracting semi-hyperbolic 1 0 node of codimension 2. PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 7 / 28 Dynamic behaviors of a modified SIR model Theorem 0.4. If dδ >βA, E is globally asymptotically stable. 0 0 Proof. If dδ >βA, it is obvious thatR < 1 andC > 0. From Theorem 0.1 and 0.3, E is 0 0 the unique attracting node of system (5). In order to prove that the disease free equilibrium E is globally and asymptotically stable, we construct the following Liapunov function: A dS dS V…S; I† ˆ ln ‡ I: …15† d A A It is easy to discover that E ; 0 attains the global minimum of the function V(S, I), so V(S, I) > 0. Along system (5), it turns out: A bAI Vj ˆ 2A dS ‡ …d‡ a‡ m…b; I††I: …16† …5† dS d…1‡ aI † Since μ(b, I) > μ for all I 0, we have 2 3 A …bA dd †I dad I 0 0 Vj  2A dS ‡  0: …17† …5† dS d…1‡ aI † The equality V…S; I† ˆ 0 if and only if at E ; 0 . By Poincare-Bendixson theorem, theorem 0.4 is obvious. Let E…S; I† be any endemic equilibrium, one can verify that its characteristic equation can be written as …18† l tr…J †l‡ det…J † ˆ 0; E E where m m 2aI bI 1 0 tr…J † ˆ d‡ bI …d‡ a‡ m…b; I†† ; 2 2 1‡ aI 1‡ aI …b‡ I† …19† det…J † ˆ …b‡ I†f …I†: E 2 …b‡ I† …1‡ aI † Obviously, the signs of the eigenvalues are determined by f …I† and tr(J ). From Fig 2, we 0 0 know that f …I † < 0; f …I † > 0, so E is a hyperbolic saddle and E is an anti-saddle. E is an 1 2 2 1 2 attracting node or focus, if tr(J ) < 0; E is a weak focus or a center, if tr(J ) = 0; E is a repelling E 2 E 2 node or focus, if tr(J ) > 0. So we obtain the following theorem. Theorem 0.5. For system (5), there are two endemic equilibria E , E whenR < 1 and 1 2 Δ < 0. Then the low endemic equilibrium E is a hyperbolic saddle, and the higher endemic equi- 0 1 librium E is an anti-saddle. WhenR > 1 there is a unique endemic equilibrium, which is an anti-saddle. Bifurcation analysis Backward bifurcation d…m m † 1 0 Theorem 0.6. WhenR ˆ 1, system (5) undergoes backward bifurcation if b < ; and sys- bd d…m m † 1 0 tem (5) undergoes forward bifurcation if b > . bd PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 8 / 28 Dynamic behaviors of a modified SIR model Proof. For convenience of the proof, we suppose that the total number of the population is N(t). System (4) becomes the following system dI b…N I R†I b ˆ dI aI m ‡…m m † I; 0 1 0 dt 1‡ aI b‡ I dR b …20† ˆ m ‡…m m † I dR; 0 1 0 dt b‡ I > dN ˆ A dN aI: dt Let V = (I, R, N) , then the disease-free equilibrium is V ˆ 0; 0; and we can write Eq (20) in vector forms as: …21† V ˆ H…V†…V V †; where 0 1 b…N I R† d a m…b; I† 0 0 B 2 C 1‡ aI B C B C H…V† ˆ : …22† B C m…b; I† d 0 B C @ A a 0 d Then, 0 1 d …R 1† 0 0 1 0 B C B C H…V † ˆ m1 d 0 : …23† B C @ A a 0 d We know that the dominant eigenvalue of H(V ) is zero, ifR ˆ 1. It is well known that we can decompose a neighborhood of the disease-free state into stable manifold W and a center manifold W . 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 is tangential to the 0 c eigenvector V at V . Thus, we can assume that W has the following form: c 0 W ˆ fV ‡ aV ‡ Z…a† : V  Z…a† ˆ 0; a  a  a g; …24† 0 0 0 where V is the dominant left eigenvector of H(V ),α > 0 is a constant, and Z: [−α ,α ]! 0 0 0 0 Ran(H(V )) satisfies: Z…0† ˆ Z…0† ˆ 0: …25† da In other words, W 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 V…t† ˆ V ‡ a…t†V ‡ Z…a…t††; …26† PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 9 / 28 Dynamic behaviors of a modified SIR model since W is an invariant, from Eq (21) we have a _…t†V ‡ Z…a…t†† ˆ V…t† dt ˆ H…V…t††‰V…t† V Š 0 0 ˆ H…V ‡ a…t†V ‡ Z…a…t†††‰a…t†V ‡ Z…a…t††Š; Multiplying both sides of the above equation by V and using the following equations: …27† V  Z…a…t†† ˆ 0; V H…V † ˆ 0; V  V ˆ 1; dt and Z(α) = O(α ) we can get that 0 0 a ˆ V  H…V ‡ aV ‡ Z…a††‰aV ‡ Z…a†Š 0 0 3 ˆ V  H…V ‡ aV †‰aV ‡ Z…a†Š‡ O…a † 0 0 3 ˆ V …H…V ‡ aV † H…V ††‰aV ‡ Z…a†Š‡ O…a †: 0 0 0 0 3 Note that [H(V +αV ) − H(V )] is of orderα, then [H(V +αV ) − H(V )]Z(α) is O(α ) and 0 0 0 0 we get that 0 0 3 a _ ˆ aV ‰H…V ‡ aV † H…V †ŠV ‡ O…a †: …28† 0 0 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: 0 2 3 a _ ˆ V  H V a ‡ O…a †; …29† where dH…V ‡ aV † @H H ˆ j ˆ V j ; …30† aˆ0 VˆV da @V which gives the rate of change of the vector field as the disease invades. Hence, the number h ˆ V H V …31† determines whether the disease can invade whenR ˆ 1, and hence gives the sign of the bifur- cation. For our system, by computation we can get the V and V as follows: m a 1 0 0 1 0 0 V ˆ I ; I ; I ; V ˆ ; 0; 0 ; d d I and 0 0 0 m1 0 a 0 0 0 H ˆ H I ‡ H I H I ; 1 2 3 d d where 0 1 m m 1 0 0 1 0 1 b‡ 0 0 b 0 0 b 0 0 B C B C B C B C 0 B C 0 B C 0 B C m m H ˆ 1 0 ; H ˆ 0 0 0 ; H ˆ 0 0 0 : …32† B C B C B C 1 2 3 0 0 B C @ A @ A @ A 0 0 0 0 0 0 0 0 0 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 10 / 28 Dynamic behaviors of a modified SIR model Then we can get that h ˆ V H V m m b…d‡ a‡ m † 1 0 1 0 ˆ I : …33† b d According to [27], we know that system (5) undergoes backward bifurcation, when h > 0, i.e., d…m m † d…m m † 1 0 1 0 b < ; and system (5) undergoes forward bifurcation, when h < 0, i.e., b > . bd1 bd1 Proposition 0.7. WhenR passes through R and tr(I )6ˆ 0, system (5) with a! 0 under- goes a saddle-node bifurcation. WhenR ˆ R , E is a saddle-node if tr(I )6ˆ 0, and E is a cusp 0 0 if tr(I ) = 0. Proof. As a! 0,Δ (a)!Δ (0).Δ (0) = 0 means thatR ˆ R and the two endemic equilib- 0 0 0 0 0 ria E and E coalesce at E . Two eigenvalues of Jacobian matrix J(E ) are 0 and tr(I ) for system 1 2 (5). 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 b …m m †…bb d† 2 3 0 1 2 < X ˆ X ‡ XO…jYj†‡ O…jYj ;jX; Yj † …b‡ I † jTj …34† Y ˆ tr…I †Y ‡ O…jX; Yj † Where T is the non-singular transformation to diagonalize the linear part. Since b < , E is a saddle-node if tr(I )6ˆ 0. Combined with Theorem 0.1, system (5) undergoes saddle-node bifurcation whenR passes through the critical valueR , 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 numberR ˆ 1 defines a straight line C in the (β, b) plane, dd C : b ˆ : C ˆ 0 also defines one branch of the hyperbolic C (see Fig 3), bA dd C : b ˆ f …b† ˆ : B C bd dd A…m m † 1 1 0 The branch of C intersects with C at the point K ; and withβÐaxis at the B 0 2 A d dd point K ; 0 . It is easily found that f is an increasing convex function ofβ in the first quad- rant. Let dd A…m m † ‡ 1 1 0 C ˆ f…b; b†jb ˆ ; b > g; 0 2 dd A…m m † 1 1 0 C ˆ f…b; b†jb ˆ ; b < g; 0 2 A d ‡ ‡ then C ˆ C [ C [ K . Here C and C are two branches of C joint at point K. 0 0 0 0 0 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 11 / 28 Dynamic behaviors of a modified SIR model Fig 3. The bifurcation curves in (β, b) for system (5) when a6ˆ 0. https://doi.org/10.1371/journal.pone.0175789.g003 Define the curveΔ (β, b) = 0 as C , one can verify that 0 D D …K † ˆ 0; D …K† ˆ 0; 0 0 @b @b j 0 ˆ 0; j ˆ 1: K K @b @b Hence, the curve C is tangent to the curve C at the point K and theβÐaxis at the point K dd dd 0 1 when < b < . A A dd If b ˆ , 2 2 3d …bd ‡ A…d d †† g…b† 0 1 D …b; b† ˆ D …b† ˆ ; 0 0 where 2 2 2 2 2 2 2 g…b† ˆ A a d b 2Aad d b d …4A a…d d † d d †: 1 0 1 0 0 1 0 1 4 3 Denote the discrimination of g(b) = 0 as D ˆ 16A a d d …d d † < 0. Hence the equation 2 0 0 1 A…m m † 1 0 Δ (b) = 0 has a unique real solution, b ˆ , which means that K is the only point at which 0 2 the curve C is tangent to the curve C . If b = 0, D …b; 0† ˆ D …b† ˆ 3d …dd bA† g …b†; 0 0 0 0 1 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 12 / 28 Dynamic behaviors of a modified SIR model where g …b† ˆ d b ‡ 4Aadb 4ad d : 1 0 0 Through computing, we find that equationΔ (β, 0) = 0 has three real solutions q 2 2 2adA 2d A a ‡ ad dd b ˆ ; b ˆ : 0 1;2 A d It is easy to verifyβ >β , so the curve C will not intersect with the abscissa axis when 0 1 dd dd 0 1 b 2 ; . A A For the curve C , note bA dd 12…ad…bA dd †=b‡ bd † …bA dd †…dd bA† 0 0 0 0 1 D b; ˆ bd bd 1 1 …35† 2 2 2 2 2 81a d d …bA dd † …dd bA† 0 0 1 ‡ : 2 2 b d dd dd bA dd 0 1 0 Obviously, if b 2 ; , then D b; > 0. Hence, the curve C is located above the 0 B A A bd dd dd 0 1 curve C for b 2 ; . 0 A A Based on the above the discussion and Theorem 0.1, if we define dd D ˆ f…b; b†jb > ; b > 0g; dd dd A…m m † 0 1 1 0 D ˆ f…b; b†j < b < ; 0 < b < ; D …b; b† < 0g; 2 2 0 A A then there is one endemic equilibria in the region D and two endemic equilibria in the region dd dd D . System (5) undergoes saddle-node bifurcation on the cure C when b 2 ; . The A A backward bifurcation occurs on the C and forward bifurcation occurs on the C . The pitch- 0 0 fork bifurcation occurs when transversally passing through the curve C at the point K. Espe- cially, 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Δ (β, b) = 0, p bA…d d †‡ d …dd bA† 2 bAd …d d †…dd bA† 1 0 0 1 0 1 0 1 …36† b ˆ f …b† ˆ : D 2 bd 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 ) = 0. Note that from Eq (19) we can rewrite tr(J ) as h …I† tr…J † ˆ ; E 2 2 …b‡ I† …1‡ aI † where 4 3 2 2 h …I† ˆ b I ‡ b I ‡ b I ‡ b I ‡ b d; …37† 1 4 3 2 1 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 13 / 28 Dynamic behaviors of a modified SIR model with b ˆ a…3d‡ 2m ‡ 2a†; 4 0 b ˆ a…6bd‡ 3bm ‡ bm ‡ 4ba†‡ b; 3 0 1 2 2 2 b ˆ a…3b d‡ 2b m ‡ 2b a†‡ 2bb‡ d; 2 1 b ˆ b…bb‡ 2d‡ m m †: 1 0 1 m m 2d 1 0 Here, h …I† is a quartic equation. Since b , b and b are non-negative, if b > , 2 3 4 m m 2d 1 0 h (I) > 0, in order to make sure that h (I) = 0 has a positive root, we must have b < . 1 1 This means sufficient number of hospital beds excludes the possibility of the disease oscilla- tion. From the expression of h …I† in Eq (37), it is not an easy task to study the Hopf bifurca- tion 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 3 2 2 h …I† ˆ bI ‡…d‡ 2bb†I ‡ b I ‡ b d: 1 1 One can verify the following lemma m m 2d 1 0 Lemma 0.8. For any positive equilibrium, if b  , we alway have tr(J )<0. In order to study Hopf bifurcation and Bogdanov-Takens bifurcation, we will assume that m m 2d 1 0 b < . Denote the discrimination of h …I† asΔ . Since b < 0, function h …I† must have 1 1 1 1 one negative real root. As shown in Fig 4, It is not difficult to verify function h …I† has two Fig 4. Graph of h(I) with different signs ofΔ when b < 0. When I = H , I = H or I = H = H , Hopf bifurcation occurs. BT bifurcation 1 1 2 m 2 M 2 M m of codimension 2 occurs when I* = H or I* = H and BT bifurcation of codimension 3 occurs when I* = H = H . m M m M https://doi.org/10.1371/journal.pone.0175789.g004 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 14 / 28 Dynamic behaviors of a modified SIR model 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 …I† is determined by the sign of theΔ . When exist, we denote the roots as H and H with H  H . m M M m @I @I 2 2 Lemma 0.9. > 0; 8b > 0 and < 0; 8b > 0. @b @b Proof. From the Eq (9) and the expression of I , direct calculation leads to @I @f @f dd I ‡ bdd 2 0 1 p ˆ = ˆ > 0; @b @b @I b D @I 1 @C 1 @C @D ˆ ‡p C 2B : @b 2B @b @b @b @I @C One can find that ˆ bd > 0. We will prove < 0 in two cases. IfR < 1, then @b 1 @b 0 @D ˆ dd bA > 0. Recall the analysis of the existence of the equilibria we know that the @b @I @I 2 2 C < 0, so the < 0. IfR > 1, then lim ˆ 0 and 0 b!‡1 @b @b " # 3=2 2 @I D @D @ D 2 0 0 0 ˆ 2D 2 0 2 8B @b @ b @ b ˆ 2b A…m m †…bA dd † > 0; 1 0 1 @I …b† so the is an increasing function of b with supremum 0 in the (0, +1), so for8b > 0 @b @I …b† < 0. @b Theorem 0.10. For system (5) with a = 0, generic Hopf bifurcation could occur if I = H , 2 m I = H or I = H = H . 2 M 2 m M Proof. We only need to verify the transversality condition. Letγ = tr(I )/2 be the real part of the two solutions of Eq (18), when a = 0. If I = H or I = H = H , we considerβ as the bifurcation parameter and fix all other 2 M 2 M m parameters. Then dg 1 @tr…I …b†; b† @I …b† @tr…I …b†; b† 2 2 2 ˆ ‡ ; db 2 @I @b @b ^ ^ bˆb 2 bˆb @tr…I …b†; b† 2h…H † 2 M ˆ < 0; @I ^ …b‡ H † 2 bˆb 3 2 2 @tr…I …b†; b† H ‡ 2bH ‡ b H 2 M M M ˆ < 0: @b ^ …b‡ H † bˆb @I dg From Lemma 0.9, we have j > 0, so < 0. bˆb @b db If I = H or I = H = H , we consider b as the bifurcation parameter and fix all the other 2 m 2 M m parameters. Then dg 1 @tr…I …b†; b†@I …b† @tr…I …b†; b† 2 2 2 ˆ ‡ ; db 2 @I @b @b ^ ^ bˆb 2 bˆb @tr…I …b†; b† 2h…H † 2 m ˆ > 0; @b ^ …b‡ H † bˆb m @tr…I …b†; b† …b H †…m m † 2 m 0 1 ˆ < 0: @b ^ …b‡ H † bˆb m PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 15 / 28 Dynamic behaviors of a modified SIR model Fig 5. Graphs of Bifurcation curve in parameters plane (β, b) and the phase trajectory for system (5). (a) Curve q(β, b) = 0. The green curve is supercritical Hopf bifurcation; The red curve is subcritical Hopf bifurcation.σ becomes 0 at the DH point. (b) Two limit cycles bifurcation from the weak focus E . https://doi.org/10.1371/journal.pone.0175789.g005 @I dg From Lemma 0.9, we have j < 0, so < 0 (one can verify that H < b). Then the proof @b bˆb dm 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 . If E is a weak focus, Hopf bifurcation can happen, otherwise system does not undergo 2 2 Hopf bifurcation. Because system (5) is analytic when a = 0, E can only be weak focus or cen- ter. We can distinguish these two types of singularities by calculating Lyapunov coefficients and verifying it through numerical simulation. Taking the resultant of f(I) and h …I† with respect to I, we can get the curve q(β, b) in the space (β, b), and plot the algebraic curve q(β, b) = 0 by fixing other parameters A, d, μ ,α and μ . Choose A = 3, d = 0.3,α = 0.5, μ = 1.5, μ = 3 and plot q(β, b) = 0 in the plane (β, b) as 0 0 1 shown in Fig 5(a). The green curve (δ < 0) represents supercritical Hopf bifurcation; the red curve corresponding toδ > 0 represents subcritical Hopf bifurcation. 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 (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 . Bogdanov-Takens bifurcation From Theorem 0.1 we know that the two equilibria E and E coalesce at the equilibria 1 2 E (S , I ) whenR ˆ R , if a = 0, where 0 0 S ˆ ; d‡ bI d…d‡ a‡ m † bA‡ bb…d‡ a‡ m † 0 1 I ˆ : 2b…d‡ a‡ m † We can find that det(I ) = 0 in Eq (18) ifR ˆ R . From Proposition 0.7, we know that E is a 0 0 saddle-node point if tr(I )6ˆ 0. If tr(I ) = 0, Eq (18) has a zero eigenvalue with multiplicity 2, PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 16 / 28 Dynamic behaviors of a modified SIR model 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 (I )6ˆ 0, then E is a Bogdanov-Takens point of codimension 2, and system (5) localized at E is topologically equivalent to X ˆ Y; …38† : 3 Y ˆ X ‡ Sign…h…I ††XY ‡ O…jX; Yj †: Proof. Changing the variables as x = S − S , y = I − I , then system (5) becomes dx > ˆ …d‡ bI †x bS y bxy; dt …39† dy …m m †bI 1 0 2 : ˆ bI x‡ y‡ Cy ‡ bxy‡ O…jx; yj †; dt …b‡ I † where …m m †b …m m †bI 1 0 1 0 C ˆ : …40† 2 3 …b‡ I † …b‡ I † tr(I ) = 0 and det(I ) = 0 imply that …m m †bI 2 1 0 b S I …d‡ bI † ˆ 0; d‡ bI ˆ ; …41† …b‡ I † so the generalized eigenvectors corresponding toλ = 0 of Jacobian matrix in system (5) at the point E are 0 0 …42† V ˆ … d bI ; bI † ; V ˆ …1; 0† : 1 2 Let T = (T ) = (V , V ), then under the non-singular linear transformation ij 2×2 1 2 ! ! x X ˆ T ; y Y where |T| = −βI < 0. System (39) becomes X ˆ Y ‡ a X ‡ bXY; …43† : 3 Y ˆ a X ‡ dbXY ‡ O…jX; Yj †; here …m m †b …m m †bI 2 1 0 1 0 2 …d‡ bI †b I b I 2 3 …b‡ I † …b‡ I † a ˆ ; bI bI …bI ‡ d†…bb d† a ˆ : b‡ I PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 17 / 28 Dynamic behaviors of a modified SIR model Using the near-identity transformation a X 12 2 …44† u ˆ X ; v ˆ Y ‡ a X ; and rewrite u, v into X, Y, we obtain X ˆ Y ‡ O…jX; Yj †; …45† : 3 Y ˆ M X ‡ M XY ‡ O…jX; Yj †: 21 22 To consider the sign of M , note that bI …bI ‡ d†…bb d† M ˆ a ˆ : 21 21 b‡ I d…m m † 1 0 For system (5), the condition of the existence of endemic equilibrium is b < , hence, bd M < 0. Then we will determine the sign of M by 21 22 M ˆ a ‡ 2a 22 22 11 2 …m m †b …m m †bI 2 1 0 1 0 …d‡ bI †b I b I 2 3 …b‡I † …b‡I † ˆ db 2 bI bh…I † ˆ : …b‡ I † If h (I )6ˆ 0, we make a change of coordinates and time and preserve the orientation by time M M a 21 21 22 X ! X; Y ! Y; t ! j jt …46† a a M 22 21 then system (5) is topologically equivalent to the normal form Eq (38). From Theorem 0.11, we know that if a = 0, endemic equilibrium E is a Bogdanov-Takens 0  0 point of codimension 2 whenΔ = 0, h(I ) = 0 and h (I )6ˆ 0. If h (I ) = 0, E may be a cusp of codimension 3. In [28], a generic unfolding with the parametersε = (ε ,ε ,ε ) of the codimension 3 cusp 1 2 3 singularity is C equivalent to X ˆ Y; …47† : 4 2 3 Y ˆ ε ‡ε Y ‡ε XY ‡ X X Y ‡ O…jX; Yj †: 1 2 3 About system (47), we can find that the system has no equilibrium ifε > 0. The planeε = 1 1 0 excluding the origin in the parameter space is saddle-node bifurcation surface. Whenε 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ε < 0. They can be visualized by drawing their trace on the half-sphere 2 2 2 2 S ˆ f…ε ;ε ;ε †jε < 0;ε ‡ε ‡ε ˆ s g; …48† 1 2 3 1 1 2 3 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 of homoclinic bifurcation, om a H of Hopf bifurcation and SN of double limit cycle bifurcation. The curve SN include two lc lc PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 18 / 28 Dynamic behaviors of a modified SIR model Fig 6. The bifurcation diagram of BT of codimension 3. (a) The parameter space and the trace of the bifurcation diagram on the S(  0); (b) The sign of the BT is positive if the coefficient of the term XY in the norm form is positive, otherwise it is negative [24]. https://doi.org/10.1371/journal.pone.0175789.g006 points H and H and the curve SN tangent to the curves H and H . The curves H and 2 om2 lc om 2 2 2 H touch the @S ˆ f…ε ;ε ;ε †jε ˆ 0;ε ‡ε ˆ s g with a first-order tangency at the om 1 2 3 1 2 3 + − + − 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 near the BT and an unique stable limit cycle between H and H near the om om BT . In the curved triangle CH H the system has two limit cycles, the inner one unstable 2 om2 and the outer one stable. These two limit cycles coalesce when theε crosses over the curve SN . On the SN there exists a unique semistable limit cycle. The more interpretation can be lc lc 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 The- orem 0.11, we make a series of changes of variables and time, so there will be different situa- tions with different signs of h (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, μ = 1.5, μ = 3, a = 0, we find that (β, b) = (0.367004, 0 1 0.183323) satisfying the conditions in Theorem 0.11, then we use (β, b) to unfold the Bogda- nov-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 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 19 / 28 Dynamic behaviors of a modified SIR model + − Fig 7. The bifurcation diagram in plane (β, 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 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 is a unstable focus. In the IV, the phase dia- gram 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 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 dia- gram 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 Fig 8(b) and 8(c), we can observe that the homoclinic loops are located in IV. Otherwise, from Fig 8(d) and 8(e), we can obtain that the homoclinic loops are located in II which is in the small neighborhood of BT . Hence, the Hopf bifurcation curve and homoclinic loops switch their positions at some point C. In order to figure out the relative positions of C, H and Hom , as shown in Fig 9 we 2 2 change the value ofβ and plot the bifurcation diagram on (β, b) with different b. PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 20 / 28 Dynamic behaviors of a modified SIR model Fig 8. The phase diagram of system (5). The blue curve represents unstable manifold, green curve represents stable manifold. (a) b = 0.1,β = 0.339; (b) b = 0.1,β = 0.340; (c) b = 0.1,β = 0.3415; (d) b = 0.18,β = 0.367; (e) b = 0.18,β = 0.3737; (f) b = 0.21,β = 0.3815; (g) b = 0.21,β = 0.4; (h) b = 0.1587,β = 0.3683. In the (b), there is an unstable limit cycle marked black curve near the epidemic equilibrium E . In the (e) and (f), there is a stable limit cycle marked red curve. In the (h), we find that there are two limit cycle, the small one is unstable, the another one is stable. https://doi.org/10.1371/journal.pone.0175789.g008 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 21 / 28 Dynamic behaviors of a modified SIR model Fig 9. Bifurcation diagram in (β, I) with different b. The red dash(solid) represents unstable epidemic equilibrium(limit cycle). The blue curve represents stable epidemic equilibrium or limit cycle. (a) b = 0.145; (b) b = 0.14. https://doi.org/10.1371/journal.pone.0175789.g009 By simulations, we get the case (a) and (b) in Fig 9, we now know that two limit cycles bifur- cating 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 > b > b H C H 2 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 bifurca- tion 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 behav- ior 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 a6ˆ 0. In the simulation, the parameters excluding a are the same as the simulation setting. From Fig 12, we find that the region D 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 A…m m † 1 0 0.3, we get the necessary and sufficient condition of backward bifurcation is b < when A…m m † dd 1 0 1 R ˆ 1, which means that we can eliminate the disease if b > and b < i.e we need m m 2d 1 0 enough number of hospital beds. From the Lemma 0.8, we know that if b > , system (5) will not have periodic solution, and the endemic equilibrium E is stable. We then discuss Hopf bifurcation and BT bifurcation for system (5) and present in details about these bifurca- tions 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 (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) ifR > 1. If R < 1, from Fig 3 and Eq (36), we know that if b > b ˆ f we can eliminate the disease, and 0 c D PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 22 / 28 Dynamic behaviors of a modified SIR model Fig 10. Bifurcation digram near the Bogdanov-Takens. https://doi.org/10.1371/journal.pone.0175789.g010 Fig 11. Phase portraits for parameters in different regions of Fig 10. https://doi.org/10.1371/journal.pone.0175789.g011 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 23 / 28 Dynamic behaviors of a modified SIR model Fig 12. The curve q(β, b) = 0 with different values of a. The blue curve, yellow curve, red curve and green curve are drawn according to a = 0, a = 0.5, a = 1 and a = 2 respectively. https://doi.org/10.1371/journal.pone.0175789.g012 these rich dynamics finally disappear through the saddle-node bifurcation when b = b as shown in Figs 7 and 13(b) and 13(c). For system (3) with a6ˆ 0, taking A = 3, d = 0.3,β = 0.5, a = 0.2, μ = 3.1728627 and μ = 1.5 1 0 we get the bifurcation diagram with different values for c as shown as Figs 14 and 15. In Fig 14, the types of BT-bifurcation are the same, however, there are two types of BT bifurcations in the Fig 15. Fig 13. Bifurcation in the plane (b, I) with differentμ .β = 0.39, a = 0. (a)R ˆ 1:02; (b)R ˆ 0:98; (c)R ˆ 0:95. 1 0 0 0 https://doi.org/10.1371/journal.pone.0175789.g013 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 24 / 28 Dynamic behaviors of a modified SIR model Fig 14. The bifurcation diagram of system (3) in parameters plane (β, b). A = 3, d = 0.3,β = 0.5, c = 0.185, a = 0.2,μ = 3.1728627,μ = 1 0 ‡ ‡ 1.5,BT …0:353073; 0:0925301†;BT …0:375; 0:137406†. 1 2 https://doi.org/10.1371/journal.pone.0175789.g014 Fig 15. The bifurcation diagram of system (3) in parameters plane (β, b). A = 3, d = 0.3,β = 0.5, c = 0.1, a = 0.2,μ = 3.1728627,μ = 1 0 + − 1.5, BT (0.337066, 0.072821), BT (0.382572, 0.172627). https://doi.org/10.1371/journal.pone.0175789.g015 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 25 / 28 Dynamic behaviors of a modified SIR model Fig 16. Phase portraits of system (3) in the plane (S, I) with different d. A = 3,β = 0.375,α = 0.5,μ = 0.5,μ = 0.7, a = 0.7, b = 0.01, c = 0 1 −1.65. (a) d = 1.49; (b) d = 1.483783; (c) d = 1.4. In case (a) (b) and (c), E is always a saddle and E is a stable node. E is a stable node in 1 0 2 case (b) and (c), while it is a unstable node in case (a). There is an unstable limit cycle near E in case (b). https://doi.org/10.1371/journal.pone.0175789.g016 In Fig 16, A = 3,β = 0.375,α = 0.5, μ = 0.5 we plot the phase portraits in plane (S, I) with bA different d, in these cases < 1, and find that there is an unstable limit cycle near the E d…d‡a‡m † when d = 1.483783. From the above stimulation, we know that Therorem 0.4 is not ture when p 2 a < c < 0. According to an early SIR model with nonmonotone incidence rate in the literature [19], the dynamics of the system are completely determined byR , which means that the disease will be eliminated ifR < 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 thatR is not enough to determine the dynamic behavior in system (5). By simulations, we predict that there would exist b in system (3)? which has the same role as b . Hopefully we can explore more 1c c relationships between the intervention actions, hospitalization conditions and spread of dis- eases, to provide the guidelines for public and desicion makers. Acknowledgments The authors would like to thank anonymous referees for very helpful suggestions and com- ments which led to improvements of our original manuscript. Author Contributions Conceptualization: GHL YXZ. Formal analysis: GHL YXZ. Writing ± original draft: GHL YXZ. Writing ± review & editing: GHL YXZ. References 1. Brauer F, Castillo-Chavez C. Mathematical models in population biology and epidemiology[M]. New York: Springer, 2001. 2. Kermack W O, McKendrick A G. A contribution to the mathematical theory of epidemics[C] Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Soci- ety, 1927, 115(772): 700±721. PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 26 / 28 Dynamic behaviors of a modified SIR model 3. Hethcote H W. The mathematics of infectious diseases[J]. SIAM review, 2000, 42(4): 599±653. https:// doi.org/10.1137/S0036144500371907 4. Gumel A B, Ruan S, Day T, et al. Modelling strategies for controlling SARS outbreaks[J]. Proceedings of the Royal Society of London B: Biological Sciences, 2004, 271(1554): 2223±2232. https://doi.org/10. 1098/rspb.2004.2800 5. Li Li. Monthly periodic outbreak of hemorrhagic fever with renal syndrome in China[J]. Journal of Biologi- cal Systems, 2016, 24, 519±533. https://doi.org/10.1142/S0218339016500261 6. Li Li. Patch invasion in a spatial epidemic model[J]. Applied Mathematics and Computation, 2015, 258, 342±349. https://doi.org/10.1016/j.amc.2015.02.006 7. Gui-Quan Sun, et al, Transmission dynamics of cholera: Mathematical modeling and control strategies [J]. Commun. Nonlinear Sci. Numer. Simulat., 2017, 45, 235±244. https://doi.org/10.1016/j.cnsns. 2016.10.007 8. Ming-Tao Li, Jin Z., Sun G., Zhang J.. Modeling direct and indirect disease transmission using multi- group model[J]. J. Math. Anal. Appl., 2017, 446, 1292±1309. https://doi.org/10.1016/j.jmaa.2016.09. 9. Ming-Tao Li, Sun G., Wu Y., Zhang J. and Jin Z.. Transmission dynamics of a multi-group brucellosis model with mixed cross infection in public farm[J]. Applied Mathematics and Computation, 2014, 237, 582±594. https://doi.org/10.1016/j.amc.2014.03.094 10. Gui-Quan Sun, Zhang Z., Global stability for a sheep brucellosis model with immigration[J], Applied Mathematics and Computation, 2014, 246, 336±345. https://doi.org/10.1016/j.amc.2014.08.028 11. Zhi-Qiang Xia, et al, Modeling the transmission dynamics of Ebola virus disease in Liberia[J]. Scientific Reports, 2015, 5, 13857. https://doi.org/10.1038/srep13857 12. Yan-Fang Wu, Li M., Sun G.. Asymptotic analysis of schistosomiasis persistence in models with general functions[J]. Journal of the Franklin Institute, 2016, 353, 4772±4784. https://doi.org/10.1016/j.jfranklin. 2016.09.012 13. Zhang J, Lou J, Ma Z, et al. A compartmental model for the analysis of SARS transmission patterns and outbreak control measures in China[J]. Applied Mathematics and Computation, 2005, 162(2): 909± 924. https://doi.org/10.1016/j.amc.2003.12.131 14. Tang S, Xiao Y, Yuan L, et al. Campus quarantine (Fengxiao) for curbing emergent infectious diseases: lessons from mitigating A/H1N1 in Xi'an, China[J]. Journal of theoretical biology, 2012, 295: 47±58. https://doi.org/10.1016/j.jtbi.2011.10.035 PMID: 22079943 15. Xiao Y, Tang S, Wu J. Media impact switching surface during an infectious disease outbreak[J]. Scien- tific reports, 2015, 5. 16. Liu R, Wu J, Zhu H. Media/psychological impact on multiple outbreaks of emerging infectious diseases [J]. Computational and Mathematical Methods in Medicine, 2007, 8(3): 153±164. https://doi.org/10. 1080/17486700701425870 17. Cui J, Sun Y, Zhu H. The impact of media on the control of infectious diseases[J]. Journal of Dynamics and Differential Equations, 2008, 20(1): 31±53. https://doi.org/10.1007/s10884-007-9075-0 18. Cui J, Tao X, Zhu H. An SIS infection model incorporating media coverage[J]. Rocky Mountain J. Math, 2008, 38(5): 1323±1334. https://doi.org/10.1216/RMJ-2008-38-5-1323 19. Xiao D, Ruan S. Global analysis of an epidemic model with nonmonotone incidence rate[J]. Mathemati- cal biosciences, 2007, 208(2): 419±429. https://doi.org/10.1016/j.mbs.2006.09.025 PMID: 17303186 20. Liu W, Levin S A, Iwasa Y. Influence of nonlinear incidence rates upon the behavior of SIRS epidemio- logical models[J]. Journal of mathematical biology, 1986, 23(2): 187±204. https://doi.org/10.1007/ BF00276956 PMID: 3958634 21. Liu W, Hethcote H W, Levin S A. Dynamical behavior of epidemiological models with nonlinear inci- dence rates[J]. Journal of mathematical biology, 1987, 25(4): 359±380. https://doi.org/10.1007/ BF00277162 PMID: 3668394 22. Wang W, Ruan S. Bifurcations in an epidemic model with constant removal rate of the infectives[J]. Journal of Mathematical Analysis and Applications, 2004, 291(2): 775±793. https://doi.org/10.1016/j. jmaa.2003.11.043 23. Wang W. Backward bifurcation of an epidemic model with treatment[J]. Mathematical biosciences, 2006, 201(1): 58±71. https://doi.org/10.1016/j.mbs.2005.12.022 PMID: 16466756 24. Shan C, Zhu H. Bifurcations and complex dynamics of an SIR model with the impact of the number of hospital beds[J]. Journal of Differential Equations, 2014, 257(5): 1662±1688. https://doi.org/10.1016/j. jde.2014.05.030 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 27 / 28 Dynamic behaviors of a modified SIR model 25. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[J]. Mathematical biosciences, 2002, 180(1): 29±48. https://doi.org/10.1016/S0025-5564(02)00108-6 PMID: 12387915 26. Huang J, Ruan S, Song J. Bifurcations in a predator?Cprey system of Leslie type with generalized Hol- ling type III functional response[J]. Journal of Differential Equations, 2014, 257(6): 1721±1752. https:// doi.org/10.1016/j.jde.2014.04.024 27. Dushoff F, Huang W, Castillo-Chavez C. Backward bifurcation and catastrophe in simple models of fatal disease[J]. Journal of mathematical biology, 1998, 36(3): 227±248. https://doi.org/10.1007/ s002850050099 PMID: 9528119 28. Dumortier F, Roussarie R, Sotomayor J. Generic 3-parameter families of vector fields on the plane, unfolding a singularity with nilpotent linear part. The cusp case of codimension 3[J]. Ergodic theory and dynamical systems, 1987, 7(03): 375±413. https://doi.org/10.1017/S0143385700004119 PLOS ONE | https://doi.org/10.1371/journal.pone.0175789 April 20, 2017 28 / 28

Journal

PLoS ONEPubmed Central

Published: Apr 20, 2017

References