http://orcid.org/0000-0001-7935-8950
Figures
Abstract
The transmission of infectious diseases has been studied by mathematical methods since 1760s, among which SIR model shows its advantage in its epidemiological description of spread mechanisms. Here we established a modified SIR model with nonlinear incidence and recovery rates, to understand the influence by any government intervention and hospitalization condition variation in the spread of diseases. By analyzing the existence and stability of the equilibria, we found that the basic reproduction number is not a threshold parameter, and our model undergoes backward bifurcation when there is limited number of hospital beds. When the saturated coefficient a is set to zero, it is discovered that the model undergoes the Saddle-Node bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation of codimension 2. The bifurcation diagram can further be drawn near the cusp type of the Bogdanov-Takens bifurcation of codimension 3 by numerical simulation. We also found a critical value of the hospital beds bc at
and sufficiently small a, which suggests that the disease can be eliminated at the hospitals where the number of beds is larger than bc. The same dynamic behaviors exist even when a ≠ 0. Therefore, it can be concluded that a sufficient number of the beds is critical to control the epidemic.
Citation: Li G-H, Zhang Y-X (2017) Dynamic behaviors of a modified SIR model in epidemic diseases using nonlinear incidence and recovery rates. PLoS ONE 12(4): e0175789. https://doi.org/10.1371/journal.pone.0175789
Editor: Gui-Quan Sun, Shanxi University, CHINA
Received: February 14, 2017; Accepted: March 31, 2017; Published: April 20, 2017
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.
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: This work is supported by the National Science Foundation of China (11331009), Research Project Supported by Fund Program for the Scientific Activities of Selected Returned Overseas Professionals in Shanxi Province and the National Sciences Foundation of Shanxi Province (2015011009).
Competing interests: The authors have declared that no competing interests exist.
Introduction
Since the development of the first dynamic model of smallpox by Bernoulli in 1760, various mathematical models have been employed to study infectious diseases [1] in order to reveal the underlying spread mechanisms that influence the transmission and control of these diseases. Among them, Kermack and Mckendrick [2] initiated a famous SIR type of compartmental model in 1927 for the plague studies in Mumbai, and succeed in unveiling its epidemiological transition. Since then, mathematical modeling has become an important tool to study the transmission and spread of epidemic diseases.
In the modeling of infectious diseases, the incidence function is one of the important factors to decide the dynamics of epidemic models. Bilinear and standard incidence rates, both monotonically increasing functions of the total of infected class, have been frequently used in early epidemic models [3]. In those models, the dynamics of models are relatively simple and almost determined by the basic reproduction number : the disease will be eliminated if
, otherwise, the disease will persist. However, intervention strategies, such as isolation, quarantine, mask-wearing and medical report about emerging infectious diseases, play an vital role in controlling the spread, sometimes contributing to the eradication of diseases. For instance, the SARS in 2003 and novel influenza pandemic in 2009 have been well controlled by taking these intervention actions [4–15]. Hence, it is essential to expand the modeling studies to the investigation of the combined effects of these major intervention strategies. The generalized models will provide further understanding of the transmission mechanisms, and modify guidelines for public health in control of the spread of infectious diseases.
In recent years, a number of compartmental models have been formulated to explore the impact of intervention strategies on the transmission dynamics of infectious diseases. If denote the total number of hospital individuals, exposed and infectious as N, E and I respectively, Liu et al [16], Cui et al [17, 18] used βe−mI, and c1 − c2 f(I) to study the impact of media coverage on the dynamics of epidemic models, respectively. However, people may adjust their behaviors according to these government intervention. Therefore, Ruan and Xiao [19] set incidence function in the form of
(a special case of
[20, 21]) to include the above “psychological” effect: when the number of infectious individuals increases and is reported through social media, the susceptible individuals will stay alert spontaneously to reduce any unnecessary contact with others, thus lowering the contact and transmission incidence.
On the other hand, medical treatments, determining how well the diseases are controlled, are normally expressed as constant recovery rates in the current models. These recovery rates depend on various health systems and hospitalization conditions, such as the capacity of the hospitals and effectiveness of the treatments. Advanced models (see [22–24]) started to corporate the limited medical resources into the spread dynamics of infectious diseases. In the literature [22], Wang and Ruan first introduced a piece-wise treatment function in an SIR model,
where the maximal treatment capacity was used to cure infectives so that the epidemic of disease can be controlled. This situation occurs only if the infectious disease needs to be eliminated due to its threats to public. They discovered that the model undergoes Saddle-Node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation, standing for the collision of two equilibria, the existance of periodic diseases, and two varying parameters in system, respectively. Wang [23] further modified the treatment rate to be proportinal to the number of infectives before the capacity of hospital was reached, by
(1)
The model was then found to perform backward bifurcation [23], indicating that the basic reproduction number was no longer a threshold.
In common hospital settings, the number of beds is an indicator of health resources, particularly the medical treatments of the infectives. Under this consideration, Shan and Zhu [24] defined the recovery rate as a function of b, the number of hospital beds, and I, the number of infectives.
where μ0 is the minimum per capita recovery rate, and μ1 the maximum per capita recovery rate. They chose the standard incidence rate and discovered the complicated dynamics including Saddle-Node bifurcation, Backward bifurcation and Bogdanov-Takens bifurcation of codimension 3, which means that the recovery rate contributes to the rich dynamics of epidemic models.
Our strategy thus becomes, both government intervention and hospitalization condition need to be incorporated to achieve a better control of the emerging infectious. Therefore, the incidence rate is expressed as
where a is positive constant and
(so that aI2 + cI + 1 > 0 for all I > 0 and hence β(I) > 0 for all I > 0). When the threshold of the number of infected individuals I* is reached, the contact transmission rate starts to decrease as the number of infected individuals grows. As shown in Fig 1, the incidence β(I) increases to its maximum and then decreases to zero as I tends to infinity, which explains the phenomenon where the rate of contacting between infected I and susceptible S decreases after government intervention. We use the same expression of hospitalization conditions as the literature [24], and the following model is then established,
(3)
where A is the recruitment rate of the susceptible population, d the natural death of the population, and α the per capita disease-induced death rate, respectively.
(a) c ≥ 0; (b) .
For system (3), the cone R3+ is a positive invariant. The C1 smoothness of the right side of system (3) implies the local existence and uniqueness of the solution with initial values in R3+. Adding up the three equations of system (3), we get N′(t) = A − dN − αI. Therefore, all solutions in the first octant approach, entering or staying inside the set, are defined by
This paper will be organized as follows. In section 2, we study the existence of the equilibria of our model. In section 3, we study the stability of the equilibria. In section 4, we examine the dynamics of the model by first looking at the backward bifurcation of system, then the much complicated Hopf bifurcation and Bogdanov-Takens bifurcation of codimension 2 and 3. We summarize our results and discuss the epidemiological significance of the number of hospital beds and intervention strategies in section 5.
Existence of equilibria
For simplicity we will focus on the case when c = 0. If c ≠ 0, but c is in small neighborhood of zero, the behaviors of model still exist. Our model thus becomes,
(4)
Since the first two equations are independent of the third, it suffices to consider the first two equations. Thus, we will focus on the reduced model in the following discussions.
We find equilibria by setting the right hand of system (5) equal to zero:
(6)
Obviously, a trivial solution of Eq (6) is
, a disease free equilibrium(DFE). For E0, by using the formula in [25], one can calculate the reproduction number
(7)
For any positive equilibrium E(S, I), also called endemic equilibrium, when exists, its S and I coordinates satisfy
(8)
where the I coordinate will be the positive root of the following cubic equation
(9)
where
Let
Denote Δ0 the discriminant of f(I) with respect to I, then
Note that
and
. As shown in Fig 2, we have the following cases about the positive roots of f(I):
- Case 1:.
. In this case,
. 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:.
. In this case,
. If
, equation f(I) = 0 has no positive solution (see Fig 1(a)). If
, 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 < 0, there are two positive solutions of the equation.
- (b). Δ0 = 0, there is unique positive solution of the equation.
- (c). Δ0 > 0, we find that Eq (9) has no positive solution.
- Case 3:.
, Eq (9) becomes
(10) If
, Eq (9) has a unique positive root. If
, Eq (9) has no positive root. Note that
means
. Thus, we get the following theorem about the equilibrium of the model.
(a) ; (b)
; (c)
; (d)
.
Theorem 0.1. For system (5) with positive parameters,
- (1). the disease-free equilibrium E0 always exists,
- (2). when
, system has a unique endemic equilibrium,
- (3). when
, and
- (a).
, system (5) does not have endemic equilibrium.
- (b).
, Δ0 < 0, there exist two endemic equilibria
and
, and when Δ0 = 0 these two endemic equilibria coalesce into the same endemic equilibrium E*; otherwise system (5) has no endemic equilibrium.
- (a).
- (4). when
, there exists a unique endemic equilibrium if and only if
i.e.
; otherwise there is no endemic equilibrium.
Remark 0.2. By calculation, we get that
(11)
When a is sufficiently small and tends to zero, the sign of Δ0 will be determined by the zero power of a. Therefore, Δ0(a) → Δ0(0) = −3β2δ0Δ as a → 0. In addition, Δ = 0 if and only if
. Here
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) = (jij)2×2, i, j = 1, 2, then a straightforward calculation gives
(12)
Firstly, we present a theorem about the disease-free equilibrium E0(A/d, 0).
Theorem 0.3. E0 is an attracting node if , and hyperbolic saddle if
. When
, E0 is a saddle-node of codimension 1 if
and attracting semi-hyperbolic node of codimension 2 if
.
Proof. For system (5), −d and are two eigenvalues of J(E0). So, E0 is an attracting node if
, and unstable if
. When
, the second eigenvalue becomes zero. In order to analyze the behavior of E0, we linearize system (5) and use the transformation of
, Y = I,
(13)
where P(X, Y) and Q(X, Y) represent the higher order terms. Obviously, E0 is a saddle-node if
. Otherwise, i.e.,
, applying the center manifold theorem, system (5) becomes
(14)
where Q1(X, Y) represents the higher order term. Thus, E0 is an attracting semi-hyperbolic node of codimension 2.
Theorem 0.4. If dδ0 > βA, E0 is globally asymptotically stable.
Proof. If dδ0 > βA, it is obvious that and
. From Theorem 0.1 and 0.3, E0 is the unique attracting node of system (5). In order to prove that the disease free equilibrium E0 is globally and asymptotically stable, we construct the following Liapunov function:
(15)
It is easy to discover that
attains the global minimum of the function V(S, I), so V(S, I) > 0. Along system (5), it turns out:
(16)
Since μ(b, I) > μ0 for all I ≥ 0, we have
(17)
The equality
if and only if at
. By Poincare-Bendixson theorem, theorem 0.4 is obvious.
Let be any endemic equilibrium, one can verify that its characteristic equation can be written as
(18)
where
(19)
Obviously, the signs of the eigenvalues are determined by
and tr(JE). From Fig 2, we know that
, so E1 is a hyperbolic saddle and E2 is an anti-saddle. E2 is an attracting node or focus, if tr(JE) < 0; E2 is a weak focus or a center, if tr(JE) = 0; E2 is a repelling node or focus, if tr(JE) > 0. So we obtain the following theorem.
Theorem 0.5. For system (5), there are two endemic equilibria E1, E2 when and Δ0 < 0. Then the low endemic equilibrium E1 is a hyperbolic saddle, and the higher endemic equilibrium E2 is an anti-saddle. When
there is a unique endemic equilibrium, which is an anti-saddle.
Bifurcation analysis
Backward bifurcation
Theorem 0.6. When , system (5) undergoes backward bifurcation if
; and system (5) undergoes forward bifurcation if
.
Proof. For convenience of the proof, we suppose that the total number of the population is N(t). System (4) becomes the following system
(20)
Let V = (I, R, N)T, then the disease-free equilibrium is
and we can write Eq (20) in vector forms as:
(21)
where
(22)
Then,
(23)
We know that the dominant eigenvalue of H(V0) is zero, if
. It is well known that we can decompose a neighborhood of the disease-free state into stable manifold Ws and a center manifold Wc. 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 Wc is tangential to the eigenvector V0 at V0. Thus, we can assume that Wc has the following form:
(24)
where V* is the dominant left eigenvector of H(V0), α0 > 0 is a constant, and Z: [−α0, α0] → Ran(H(V0)) satisfies:
(25)
In other words, Wc 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
(26)
since Wc is an invariant, from Eq (21) we have
Multiplying both sides of the above equation by V* and using the following equations:
(27)
and Z(α) = O(α2) we can get that
Note that [H(V0 + αV0) − H(V0)] is of order α, then [H(V0 + αV0) − H(V0)]Z(α) is O(α3) and we get that
(28)
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:
(29)
where
(30)
which gives the rate of change of the vector field as the disease invades. Hence, the number
(31)
determines whether the disease can invade when
, and hence gives the sign of the bifurcation. For our system, by computation we can get the V* and V0 as follows:
and
where
(32)
Then we can get that
(33)
According to [27], we know that system (5) undergoes backward bifurcation, when h > 0, i.e.,
; and system (5) undergoes forward bifurcation, when h < 0, i.e.,
.
Proposition 0.7. When passes through
and tr(I*) ≠ 0, system (5) with a → 0 undergoes a saddle-node bifurcation. When
, E* is a saddle-node if tr(I*) ≠ 0, and E* is a cusp if tr(I*) = 0.
Proof. As a → 0, Δ0(a) → Δ0(0). Δ0(0) = 0 means that and the two endemic equilibria E1 and E2 coalesce at E*. Two eigenvalues of Jacobian matrix J(E*) are 0 and tr(I*) for system (5).
If tr(I*) ≠ 0, we can linearize system (5) at the E* and diagonalize the linear part. Then we can get the following form
(34)
Where T is the non-singular transformation to diagonalize the linear part. Since
, E* is a saddle-node if tr(I*) ≠ 0. Combined with Theorem 0.1, system (5) undergoes saddle-node bifurcation when
passes through the critical value
, as a → 0. If tr(I*) = 0, E* is a cusp and we will prove it in the next section.
Based on the above analysis, we know that system (5) undergoes some bifurcation. In order to consider the impact of hospital bed number and the incidence rate on the dynamics of the model, we will chose b and β as bifurcation parameters to describe these bifurcations. The basic production number defines a straight line C0 in the (β, b) plane,
also defines one branch of the hyperbolic CB (see Fig 3),
The branch of CB intersects with C0 at the point and with β—axis at the point
. It is easily found that
is an increasing convex function of β in the first quadrant. Let
then
. Here
and
are two branches of C0 joint at point K.
Define the curve Δ0(β, b) = 0 as , one can verify that
Hence, the curve
is tangent to the curve C0 at the point K and the β—axis at the point K′ when
.
If ,
where
Denote the discrimination of g(b) = 0 as
. Hence the equation Δ0(b) = 0 has a unique real solution,
, which means that K is the only point at which the curve
is tangent to the curve C0.
If b = 0,
where
Through computing, we find that equation Δ0(β, 0) = 0 has three real solutions
It is easy to verify β0 > β1, so the curve
will not intersect with the abscissa axis when
.
For the curve CB, note
(35)
Obviously, if
, then
. Hence, the curve CB is located above the curve
for
.
Based on the above the discussion and Theorem 0.1, if we define
then there is one endemic equilibria in the region D1 and two endemic equilibria in the region D2. System (5) undergoes saddle-node bifurcation on the cure
when
. The backward bifurcation occurs on the
and forward bifurcation occurs on the
. The pitchfork bifurcation occurs when transversally passing through the curve C0 at the point K. Especially, if a = 0, system (5) has a semi-hyperbolic node of codimension 2 at the point K and we can solve b in term of β from Δ0(β, b) = 0,
(36)
Now, we discuss the Hopf bifurcation of system (5). It follows from Eq (18) that if Hopf bifurcation occurs at one endimic equilibrium , we have tr(JE) = 0. Note that from Eq (19) we can rewrite tr(JE) as
where
(37)
with
Here, is a quartic equation. Since b2, b3 and b4 are non-negative, if
, h1(I) > 0, in order to make sure that h1(I) = 0 has a positive root, we must have
. This means sufficient number of hospital beds excludes the possibility of the disease oscillation. From the expression of
in Eq (37), it is not an easy task to study the Hopf bifurcation from the polynomial Eq (37), we will study a simple case when a = 0, and give the simulations to explore the case when a is small.
When a = 0, the polynomial Eq (37) is reduced to
One can verify the following lemma
Lemma 0.8. For any positive equilibrium, if , we alway have tr(JE)<0.
In order to study Hopf bifurcation and Bogdanov-Takens bifurcation, we will assume that . Denote the discrimination of
as Δ1. Since b1 < 0, function
must have one negative real root. As shown in Fig 4, It is not difficult to verify function
has two humps which locate on the different sides of vertical axis, and the maximum is obtained on the left hump, while the minimum is obtained on the right hump.
When I2 = Hm, I2 = HM or I2 = HM = Hm, Hopf bifurcation occurs. BT bifurcation of codimension 2 occurs when I* = Hm or I* = HM and BT bifurcation of codimension 3 occurs when I* = Hm = HM.
The number of roots of function is determined by the sign of the Δ1. When exist, we denote the roots as Hm and HM with HM ≥ Hm.
Lemma 0.9. .
Proof. From the Eq (9) and the expression of I2, direct calculation leads to
One can find that
. We will prove
in two cases. If
, then
. Recall the analysis of the existence of the equilibria we know that the
, so the
. If
, then
and
so the
is an increasing function of b with supremum 0 in the (0, +∞), so for ∀b > 0
.
Theorem 0.10. For system (5) with a = 0, generic Hopf bifurcation could occur if I2 = Hm, I2 = HM or I2 = Hm = HM.
Proof. We only need to verify the transversality condition. Let γ = tr(I2)/2 be the real part of the two solutions of Eq (18), when a = 0.
If I2 = HM or I2 = HM = Hm, we consider β as the bifurcation parameter and fix all other parameters. Then
From Lemma 0.9, we have
, so
.
If I2 = Hm or I2 = HM = Hm, we consider b as the bifurcation parameter and fix all the other parameters. Then
From Lemma 0.9, we have
, so
(one can verify that Hm < b). Then the proof of theorem is completed.
The reason why we choose different parameters to unfold Hopf bifurcation in Theorem 0.10 is that the transversality condition may fail at some point if we focus on one parameter.
In order to verify that Hopf bifurcation occurs in the system, we need to know the type of E2. If E2 is a weak focus, Hopf bifurcation can happen, otherwise system does not undergo Hopf bifurcation. Because system (5) is analytic when a = 0, E2 can only be weak focus or center. We can distinguish these two types of singularities by calculating Lyapunov coefficients and verifying it through numerical simulation.
Taking the resultant of f(I) and 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, μ1, α and μ0. Choose A = 3, d = 0.3, α = 0.5, μ0 = 1.5, μ1 = 3 and plot q(β, b) = 0 in the plane (β, b) as shown in Fig 5(a). The green curve (δ1 < 0) represents supercritical Hopf bifurcation; the red curve corresponding to δ1 > 0 represents subcritical Hopf bifurcation.
(a) Curve q(β, b) = 0. The green curve is supercritical Hopf bifurcation; The red curve is subcritical Hopf bifurcation. σ1 becomes 0 at the DH point. (b) Two limit cycles bifurcation from the weak focus E2.
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 → +∞, the trajectory starting at (9, 0.1) spirals outward to the stable limit cycle (red curve) and E2(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 E1.
Bogdanov-Takens bifurcation
From Theorem 0.1 we know that the two equilibria E1 and E2 coalesce at the equilibria E*(S*, I*) when , if a = 0, where
We can find that det(I*) = 0 in Eq (18) if
. From Proposition 0.7, we know that E* is a saddle-node point if tr(I*) ≠ 0. If tr(I*) = 0, Eq (18) has a zero eigenvalue with multiplicity 2, which suggests that system (5) may admit a Bogdanov-Takens bifurcation. Then, we give the following theorem.
Theorem 0.11. For system (5) with a = 0, suppose that Δ = 0, h(I*) = 0 and h′(I*) ≠ 0, then E* is a Bogdanov-Takens point of codimension 2, and system (5) localized at E* is topologically equivalent to
(38) Proof. Changing the variables as x = S − S*, y = I − I*, then system (5) becomes
(39)
where
(40) tr(I*) = 0 and det(I*) = 0 imply that
(41)
so the generalized eigenvectors corresponding to λ = 0 of Jacobian matrix in system (5) at the point E* are
(42)
Let T = (Tij)2×2 = (V1, V2), then under the non-singular linear transformation
where |T| = −βI* < 0. System (39) becomes
(43)
here
Using the near-identity transformation
(44)
and rewrite u, v into X, Y, we obtain
(45)
To consider the sign of M21, note that
For system (5), the condition of the existence of endemic equilibrium is
, hence, M21 < 0. Then we will determine the sign of M22 by
If h′(I*) ≠ 0, we make a change of coordinates and time and preserve the orientation by time
(46)
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 point of codimension 2 when Δ = 0, h(I*) = 0 and h′(I*) ≠ 0. If h′(I*) = 0, E* may be a cusp of codimension 3.
In [28], a generic unfolding with the parameters ε = (ε1, ε2, ε3) of the codimension 3 cusp singularity is C∞ equivalent to
(47)
About system (47), we can find that the system has no equilibrium if ε1 > 0. The plane ε1 = 0 excluding the origin in the parameter space is saddle-node bifurcation surface. When ε1 decrease from this surface, the saddle-node point of Eq (47) becomes a saddle and a node. Then the other bifurcation surfaces are situated in the half space ε1 < 0. They can be visualized by drawing their trace on the half-sphere
(48)
when σ > 0 sufficiently small (see Fig 6(a)). We recall that the bifurcation set is a ‘cone’ based on its trace with S.
(a) The parameter space and the trace of the bifurcation diagram on the S(ϵ1 ≤ 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].
In Fig 6(b), trace on the S which consists of 3 curves: a curve Hom of homoclinic bifurcation, a H of Hopf bifurcation and SNlc of double limit cycle bifurcation. The curve SNlc include two points H2 and Hom2 and the curve SNlc tangent to the curves H and Hom. The curves H and Hom touch the with a first-order tangency at the points BT+ and BT−. In the neighbourhood of the BT+ and BT−, one can find the unfolding of the cusp-singularity of codimension 2. For system (47), there exists an unstable limit cycle between H and Hom near the BT+ and an unique stable limit cycle between H and Hom near the BT−. In the curved triangle CH2Hom2 the system has two limit cycles, the inner one unstable and the outer one stable. These two limit cycles coalesce when the ε crosses over the curve SNlc. On the SNlc there exists a unique semistable limit cycle. The more interpretation can be found in literature [24, 28].
Bifurcation diagram and simulation
According to analysis and Theorem 0.11, we know that system (5) undergoes Bogdanov-Takens bifurcation of codimension 2. In this section we will choose the parameters β and b as bifurcation parameters to present the bifurcation diagram by simulations. In the proof of Theorem 0.11, we make a series of changes of variables and time, so there will be different situations with different signs of h′(I). According to the positive and negative coefficients of XY term in the normal form Eq (47), we denote the Bogdanov-Takens bifurcation of codimension 2 as BT+ and BT− respectively.
Taking A = 3, d = 0.3, α = 0.5, μ0 = 1.5, μ1 = 3, a = 0, we find that (β, b) = (0.367004, 0.183323) satisfying the conditions in Theorem 0.11, then we use (β, b) to unfold the Bogdanov-Takens bifurcation of codimension 2. By simulation, we obtain the bifurcation diagram in plane (β, b) shown as Fig 7, the blue dash (solid) curve represents saddle-node bifurcation (neutral saddle), the green (blue solid) curve represents supercritical (subcritical) Hopf 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.
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).
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 E2. 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.
In Fig 8(a), β, b are located in the area between saddle-node bifurcation and subcritical Hopf bifurcation curve and the epidemic equilibrium E2 is a unstable focus. In the IV, the phase diagram of system is one of the cases shown as (b), (c) and (h). There is an unstable limit cycle (black curve) near the epidemic equilibrium E2 in Fig 8(b) and two limit cycles in Fig 8(h) with the inner one unstable and the other one stable. When β and b are located in II, the phase diagram of system is one of the cases as shown in (d) and (e) and there is a stable limit cycle in Fig 8(e). When β and b are located in I or III, the phase portraits are similar to the cases of (f) and (g), respectively. In the case (f), system (5) has a unique epidemic equilibrium and a stable limit cycle.
In the small neighborhood of BT+, we know that the unstable limit cycle bifurcating from Hopf bifurcation curve disappears from the homoclinic loop, and from 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, H2 and Hom2, as shown in Fig 9 we change the value of β and plot the bifurcation diagram on (β, b) 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.
By simulations, we get the case (a) and (b) in Fig 9, we now know that two limit cycles bifurcating from the semi-stable cycle with one disappearing from the Homoclinic loop and the another disappearing from the Hopf bifurcation curve. We therefore obtain by the order these two limit cycle vanishing with different values of b.
Then we obtain the bifurcation diagrams of system (5) near the Bogdanov-Takens bifurcation and the phase portraits in some regions of parameters shown as Figs 10 and 11 respectively.
From the above dynamical analysis, we know that system (5) has complex dynamic behavior even though a = 0. For system (5), we also find the same phenomenon by the simulation as shown in the Fig 12 for the case a ≠ 0. In the simulation, the parameters excluding a are the same as the simulation setting. From Fig 12, we find that the region D2 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.
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.
Discussion and application
In this paper we consider the SIR model with the nonmonotone incidence rate due to the intervention strategies and nonlinear recovery rate considering the hospitalization conditions.
From Theorem 0.1, we know that system (4) undergoes backward bifurcation. In Theorem 0.3, we get the necessary and sufficient condition of backward bifurcation is when
, which means that we can eliminate the disease if
and
i.e we need enough number of hospital beds. From the Lemma 0.8, we know that if
, system (5) will not have periodic solution, and the endemic equilibrium E2 is stable. We then discuss Hopf bifurcation and BT bifurcation for system (5) and present in details about these bifurcations in the case a = 0 and present the bifurcation diagrams in Figs 8 and 10.
From the discussion we get Lemma 0.9, which implies that I2(b) is a monotone decreasing function of b. Hence, increasing the number of beds can only reduce the number of the total infected individuals, but can not eliminate the diseases as shown in Fig 13(a) if . If
, from Fig 3 and Eq (36), we know that if
we can eliminate the disease, and these rich dynamics finally disappear through the saddle-node bifurcation when b = bc as shown in Figs 7 and 13(b) and 13(c).
β = 0.39, a = 0. (a) ; (b)
; (c)
.
For system (3) with a ≠ 0, taking A = 3, d = 0.3, β = 0.5, a = 0.2, μ1 = 3.1728627 and μ0 = 1.5 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.
A = 3, d = 0.3, β = 0.5, c = 0.185, a = 0.2, μ1 = 3.1728627, μ0 = 1.5,
.
A = 3, d = 0.3, β = 0.5, c = 0.1, a = 0.2, μ1 = 3.1728627, μ0 = 1.5, BT+(0.337066, 0.072821), BT−(0.382572, 0.172627).
In Fig 16, A = 3, β = 0.375, α = 0.5, μ0 = 0.5 we plot the phase portraits in plane (S, I) with different d, in these cases , and find that there is an unstable limit cycle near the E2 when d = 1.483783. From the above stimulation, we know that Therorem 0.4 is not ture when
.
A = 3, β = 0.375, α = 0.5, μ0 = 0.5, μ1 = 0.7, a = 0.7, b = 0.01, c = −1.65. (a) d = 1.49; (b) d = 1.483783; (c) d = 1.4. In case (a) (b) and (c), E1 is always a saddle and E0 is a stable node. E2 is a stable node in case (b) and (c), while it is a unstable node in case (a). There is an unstable limit cycle near E2 in case (b).
According to an early SIR model with nonmonotone incidence rate in the literature [19], the dynamics of the system are completely determined by , which means that the disease will be eliminated if
, otherwise the disease persist. Contrasting to their work and the other results for classic epidemic models, we find that the nonlinear recovery rate is also an important factor which leads to very complicated dynamics. Moreover, we find that
is not enough to determine the dynamic behavior in system (5). By simulations, we predict that there would exist b1c in system (3)? which has the same role as bc. Hopefully we can explore more relationships between the intervention actions, hospitalization conditions and spread of diseases, to provide the guidelines for public and desicion makers.
Acknowledgments
The authors would like to thank anonymous referees for very helpful suggestions and comments 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 Society, 1927, 115(772): 700–721.
- 3. Hethcote H W. The mathematics of infectious diseases[J]. SIAM review, 2000, 42(4): 599–653.
- 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.
- 5. Li Li. Monthly periodic outbreak of hemorrhagic fever with renal syndrome in China[J]. Journal of Biological Systems, 2016, 24, 519–533.
- 6. Li Li. Patch invasion in a spatial epidemic model[J]. Applied Mathematics and Computation, 2015, 258, 342–349.
- 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.
- 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.
- 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.
- 10. Gui-Quan Sun, Zhang Z., Global stability for a sheep brucellosis model with immigration[J], Applied Mathematics and Computation, 2014, 246, 336–345.
- 11. Zhi-Qiang Xia, et al, Modeling the transmission dynamics of Ebola virus disease in Liberia[J]. Scientific Reports, 2015, 5, 13857.
- 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.
- 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.
- 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. pmid:22079943
- 15. Xiao Y, Tang S, Wu J. Media impact switching surface during an infectious disease outbreak[J]. Scientific 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.
- 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.
- 18. Cui J, Tao X, Zhu H. An SIS infection model incorporating media coverage[J]. Rocky Mountain J. Math, 2008, 38(5): 1323–1334.
- 19. Xiao D, Ruan S. Global analysis of an epidemic model with nonmonotone incidence rate[J]. Mathematical biosciences, 2007, 208(2): 419–429. pmid:17303186
- 20. Liu W, Levin S A, Iwasa Y. Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models[J]. Journal of mathematical biology, 1986, 23(2): 187–204. pmid:3958634
- 21. Liu W, Hethcote H W, Levin S A. Dynamical behavior of epidemiological models with nonlinear incidence rates[J]. Journal of mathematical biology, 1987, 25(4): 359–380. 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.
- 23. Wang W. Backward bifurcation of an epidemic model with treatment[J]. Mathematical biosciences, 2006, 201(1): 58–71. 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.
- 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. pmid:12387915
- 26. Huang J, Ruan S, Song J. Bifurcations in a predator?Cprey system of Leslie type with generalized Holling type III functional response[J]. Journal of Differential Equations, 2014, 257(6): 1721–1752.
- 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. 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.