Bifurcations of an epidemic model with non-monotonic incidence rate of saturated mass action.

In this paper, we study the bifurcations of an epidemic model with non-monotonic incidence rate of saturated mass action, which describes the psychological effects of the community on certain serious diseases when the number of infective is getting larger. By carrying out the bifurcation analysis of the model, we show that there exist some values of the model parameters such that numerous kinds of bifurcation occur for the model, such as Hopf bifurcation, Bogdanov-Takens bifurcation.

Introduction

In the mathematical modeling of disease transmission, there is a classic model proposed by Kermack and McKendrick [8] in 1927. They divided the population being studied in time t into three classes labeled S(t), I(t) and R(t), where S(t) is the number of susceptible individuals, I(t) is the number of infective individuals, and R(t) is the number of removed individuals at time t, respectively. And they assumed that a rate of contacts by an infective with a susceptible is proportional to population size with constant of proportionality. It is clear that the assumption is too simple. Later Capasso and Serio [2] introduced a saturated incidence rate g(I)S into the epidemic model after studying the cholera epidemic spread in Bari in 1973, which describes the contact between infective individuals and susceptible individuals, where kI measures the infection force of the disease and 1/(1 +  αI) measures the inhibition effect from the behavioral change of the susceptible individuals when their number increases. g(I) tends to a saturation level when I gets large.

The general incidence ratewas proposed by Liu et al. [13], [14] and used by a number of authors, see, for example, [3], [5], [7], [1], etc. For a specific nonlinear incident rateRuan and Wang in [16] studied an epidemic model with this specific nonlinear incident rate and obtained lots of interesting dynamical behaviors of the model such as a limit cycle, two limit cycles and homoclinic loop, etc. Note that these functions g(I) of the incident rate g(I)S is monotone, which means that the contact rate between infective individuals and susceptible individuals is bigger and bigger as the number of infective individuals is getting larger. However, it is not true in reality. For example, in the recent outbreak of epidemic of severe acute respiratory syndrome (SARS) the “psychological” effects on the general public (see [4], [9]), aggressive measures and policies, such as border screening, mask wearing, quarantine, isolation, etc. may tend to reduce the number of contacts per unit time as the number of infective’s were getting relatively larger. To model this phenomenon, Xiao and Ruan in [18] considered a special nonlinear incident ratewhere g(I) is non-monotone. When I is small, g(I) is increasing. Because of small number of infective’s, people may ignore the epidemic, which may lead to lots of effective contacts between the infective individuals and suspectable. And when I is getting larger and larger, g(I) is decreasing since many protection measures could be taken by the susceptible individuals. In [18], Xiao and Ruan presented the global analysis of an epidemic model with the non-monotonic incident rate and obtained either the number of infective individuals tends to zero as time evolves or the disease persists. Hence, the model cannot undergo any bifurcations.

In this paper we consider a general non-monotonic incidence ratein an epidemic model, which is described by the following form:where S(t), I(t) and R(t) denote the numbers of susceptible, infective, and recovered individuals at time t, respectively, b is the recruitment rate of the population, d is the natural death rate of the population, k is the proportionality constant, μ is the natural recovery rate of the infective individuals, δ is the rate at which recovered individuals lose immunity and return to the susceptible class, α is a positive parameter, β is a parameter such that 1 +  βI  +  αI 2  > 0 for all I  ⩾ 0, hence, .

From the standpoint of biology, we are interested only in the dynamics of system (1.1) in the first octant of . Before going to details of dynamics for (1.1), we first present the following lemma which is easily proved.

The plane S   +   I   +   R   =   b/d is an invariant manifold of system (1.1) , which is attracting in the first octant of .

This lemma implies that the limit set of system (1.1) in the first octant of locates on the plane S  +  I  +  R  =  b/d. Therefore, the dynamics of system (1.1) in the first octant of is equivalent to the following system:in the first quadrant of . System (1.1) has the disease-free equilibrium and the endemic equilibria if and only if system (1.2) has the equilibrium (0, 0) and the positive equilibria, respectively. It is clear that system (1.1) always has a disease-free equilibrium E 0  = (b/d, 0, 0) for all values of parameters. To find the endemic equilibria, we investigate the existence of the positive equilibria of system (1.2). For simplicity, we re-scale (1.2) byThen we obtainwhereIt can been seen that A, p, q, n are positive parameters, .

In this paper, we focus on studying the existence of non-hyperbolic positive equilibria of (1.3) and their bifurcations. To find the positive equilibria of system (1.3), we setwhich yieldsNote thatwhere , .

Setting R 0  =  d 1/d, we have p  −  A  ⩾ 0 if and only if R 0  ⩽ 1. LetThen we obtain the following lemma.

System (1.3) has a unique positive equilibrium E ∗ (x ∗ ,   y ∗ ) if and only if one of the following conditions holds:

R  >  1; in this case, , y ∗   =   qx ∗ .

R   =   1 and 1   +   q   +   pm  <  0; in this case, , y ∗   =   qx ∗ .

Δ   =   0 and 1   +   q   +   pm  <  0; in this case, , y ∗   =   qx ∗ .

System (1.3) has two positive equilibria E (x ,   y ) and E (x ,   y ) if and only if Δ  >  0, R  <  1 and 1   +   q   +   pm  <  0. In this case, , y   =   qx ; , y   =   qx .

The Jacobian matrix of system (1.3) at equilibrium (x,  y) isTherefore, the determinant of the matrix M 1 isIts sign is determined byAnd the trace of the matrix M 1 isthe sign of which is determined byNote that pnx 2  + (1 +  q  +  pm)x  +  p  −  A  = 0. Then we havewhereFrom qualitative analysis, we obtain the following conclusion.

The unique positive equilibrium E ∗ (x ∗ ,   y ∗ ) of system (1.3) is a degenerate equilibrium if Δ   =   0 and 1   +   q   +   pm  <  0, where , y ∗   =   qx ∗ .

The unique positive equilibrium E ∗ (x ∗ ,   y ∗ ) of system (1.3) is a center-type equilibrium if R  >  1, S (x ∗ )   =   0, where

System (1.3) has two positive equilibria E (x ,   y ) and E (x ,   y ) if Δ  >  0, R  <  1 and 1   +   q   +   pm  <  0. And further when S (x )   =   0, E (x ,   y ) is a hyperbolic saddle, E (x ,   y ) is a center-type equilibrium, where

We will show that there exist some parameters values such that (1.3) has a cusp of codimension 3, and other parameters values such that (1.3) has a multiple focus. Choosing the original parameters of the model as bifurcation parameters, we discuss what bifurcations system (1.3) can undergo.

This paper is organized as follows. In Section 2, we show that there exist some values of parameters such that the model (1.3) has a unique positive equilibrium, which is a cusp of codimension 3. And there exists a set of values of parameters such that the model undergoes Bogdanov–Takens bifurcation of codimension 2 when two parameters vary in the small neighborhood of the set of parameter values. In Section 3, we show that there exist some values of parameters such that the model (1.3) has a positive equilibrium, which is a multiple focus of codimension 1 in two cases. Choosing one parameter of the model as a bifurcation parameter, we discuss the Hopf bifurcation of the model. The paper ends with a brief conclusion.

Bogdanov–Takens bifurcation

The purpose of this section is to study if there exist some values of model parameters such that model (1.3) undergoes the Bogdanov–Takens bifurcation. From Theorem 1.3, we know that the unique equilibrium E ∗(x ∗,  y ∗) is degenerate if and only if (H1) Δ  = 0 and (H2) 1 +  q  +  pm  < 0. In order to guarantee the existence of Bogdanov–Takens bifurcation, we further assume that (H3) tr(M 1(x ∗)) = 0. Now we choose some values of parameters n, p, q, m and A such that (H1)–(H3) hold. We first show that there exists the value of model parameters such that E ∗(x ∗,  y ∗) is a cusp of codimension 3. Taking m  = −7, p  = 2 and q  = 5, we obtain that n  = 20, A  = 8/5, x ∗  = 1/10, y ∗  = 1/2 by (H1)–(H3).

When (m,   p,   q,   n,   A)   =   (−7,   2,   5,   20,   8/5), system (1.3) has a unique positive equilibrium (x ∗ ,   y ∗ )   =   (1/10,   1/2), which is a cusp of codimension 3.

Translating the unique positive equilibrium (1/10, 1/2) into the origin, we set X  =  x  − 1/10, Y  =  y  − 1/2, for simplicity, rewriting X, Y as x, y, respectively, then (1.3) becomesUsing Taylor expansion to (2.1), we obtainLet X  =  x, Y  =  x  −  y/5, and rename X, Y as x, y. Then (2.2) is transformed into

In order to obtain the canonical normal forms, we set X  =  x  − 8x 2  + 200x 3/3 − 500x 4, Y  =  y  − 8x 2, and rewriting X, Y as x, y, respectively. Then (2.3) becomesFrom theorem in [19], the equilibrium (0, 0) of (2.4) is a cusp of codimension 3, which implies the conclusion. □

From Lemma 2.1, we can see when model parameters (m,  p,  q,  n,  A) vary in the small neighborhood of (−7, 2, 5, 20, 8/5), system (1.3) may undergo the cusp bifurcation of codimension 3 and some complicated dynamics will occur. However, this analysis of the cusp bifurcation of codimension 3 is too complicated to be done. We have to discuss if there exist model parameters such that system (1.3) can undergo the cusp bifurcation of codimension 2.

Taking m  = −4, p  = 2, q  = 3, from (H1) and (H3) we have A  = 5/3, n  = 6. In this case, system (1.3) has a unique positive equilibrium (1/6, 1/2). It is easy to check that the equilibrium is a cusp of codimension 2.

In the following we will find the universal unfolding of E(1/6, 1/2) by choosing parameters p and q as bifurcation parameters in a small neighborhood of (p,  A,  q,  m,  n) = (2, 5/3, 3, −4, 6). LetAnd rewriting τ as t, then (1.3) becomeswhere ∣λ 1∣ ≪ 1 and ∣λ 2∣ ≪ 1.

If λ 1  = 0, λ 2  = 0, then (1/6, 1/2) is a cusp of codimension 2 for (2.5). Let X  =  x  − 1/6, Y  =  y  − 1/2. By the Taylor expansion, we have (for simplicity, we still use x, y as X, Y, respectively)where λ  = (λ 1,  λ 2), w 1(x,  y,  λ) is a smooth function of x, y and λ at least of order three in x and y.

Set X  =  x, Y  = −λ 1/6 + (1 + 2λ 1/3)x  −  y/3, and rewriting X, Y as x, y, respectively, we getwhere w 2(x,  y,  λ), w 3(x,  y,  λ) is a smooth function of x, y and λ at least of order three in x and y, and

By setting X  =  x  − 5x 2, Y  =  y  + (−4 − 20λ 1/3)x 2  +  w 2(x,  y,  λ), and rename X, Y as x, y, respectively, we havewhere w 4(x,  y,  λ) is a smooth function of x, y and λ at least of order three in x and y, and

Obviously, b 1 b 2  ≠ 0 if ∣λ 1∣ and ∣λ 2∣ are vary small. Thus we set , Y  =  y and rename X, Y as x, y, respectively, then we obtainwhere w 5(x,  y,  λ) is a smooth function of x, y and λ at least of order three in x and y, and

Now we make the final change of variables by , , τ  = −b 1 t/b 2. And rewriting X, Y, τ as x, y, t, respectively. We getwhere w 6(x,  y,  λ) is a smooth function of x, y and λ at least of order three in x and y, and

LetAnd after simple calculation we can obtain thatThus, by the Bogdanov and Takens bifurcation theorems, we obtain the following conclusion in a small neighborhood of (λ 1,  λ 2) = (0, 0).

System (2.5) undergoes a saddle-node bifurcation, and the saddle-node bifurcation curve is .

System (2.5) undergoes a Hopf bifurcation, and there exist some values of λ and λ such that system (2.5) has a unique stable limit cycle.

System (2.5) undergoes a Homoclinic bifurcation, and there exist some other values of λ and λ such that system (2.5) has a unique stable homoclinic loop.

For example, taking λ 1  = 0.05, λ 2  = −0.21, system (2.5) has two positive equilibria E 1 and E 2, and a stable limit cycle (SLC) (see Fig. 2.1 ), where E 1 is a saddle and E 2 is an unstable focus.

Fig. 2.1

The phase portrait of system (2.5) when λ1 = 0.05, λ2 = −0.21.

Hopf bifurcation

In this section, we will study the Hopf bifurcation of system (1.3) for some values of model parameters in two cases: (1) R 0  > 1; (2) Δ  > 0, R 0  < 1 and 1 +  q  +  pm  < 0.

In the case that R0 > 1

From Theorem 1.3, we know that when R 0  > 1, there is always an unique positive equilibrium (x ∗,  y ∗) of (1.3), and (x ∗,  y ∗) is a center when S 2(x ∗) = 0, whereNow, we setand substitute them into tr(M 1(x ∗)) = 0, thus, we get the equationwhich has two real roots A  = 12 and A  = 3738/361. Obviously,andThe two sets of parameter values satisfy the conditionsrespectively. We first discuss the local property of the unique positive equilibrium (x ∗,  y ∗) when (m,  n,  p,  q,  A) = (−3, 10, 9, 47, 12). In this case, we have x ∗  = 1/10, y ∗  = 47/10. We consider

(1/10,   47/10) is a stable weak focus of multiple one for system (3.1) .

In order to translate (x ∗,  y ∗) to the origin, we set X  =  x  − 1/10, Y  =  y  − 47/10. And we rename X, Y as x, y, respectively. Then (3.1) becomesNow, we only need to prove that the origin is a stable weak focus of multiple one of (3.2). Using Taylor expansion for (3.2), then we haveMaking the change of variables u  =  x, , , and rewriting τ as t, we obtainThus, we can get the first Liapounov constant . Therefore the origin is a stable weak focus of order one for (3.4). The conclusion follows. □

In the following, we choose A as a bifurcation parameter. LetThen we may write (3.1) as followsThe positive equilibrium of system (3.5) is and the linearizing matrix of (3.5) at E ∗ iswhereTherefore, the characteristic equation of which isand

Obviously, (i) Re  λ(ϵ 1)∣(ϵ 1  = 0) = 0, ; (ii) (Re  λ(ϵ 1))′(ϵ 1  = 0) = −15/104 ≠ 0; (iii) from Lemma 3.1, we have W 1  < 0.

Therefore, by Hopf bifurcation theory, we obtain the following result.

There exists a σ  >  0, and a function ϵ   =   ϵ (x ) defined on , which satisfies . And when ϵ   =   ϵ (x )  <  0, system (3.5) has a unique stable limit cycle which passes through (x ,   47/10).

Next we consider system (1.3) when (m,  n,  p,  q,  A) = (−3, 10, 9, 47, 3738/361):In this case system (1.3) has a unique positive equilibrium (x ∗,  y ∗), where x ∗  = 1/19, y ∗  = 47/19. By the similar arguments for the case (m,  n,  p,  q,  A) = (−3, 10, 9, 47, 12), we obtain that the first Liapounov constant of (1/19, 47/19) is . Therefore, we get

(1/19,   47/19) is a stable weak focus of order one for system (3.6) .

In the following, we still choose A as a bifurcation parameter. LetThen we can write (3.6) as follows:The positive equilibrium of system (3.7) is and the linearizing matrix of (3.7) at E ∗ iswhere

Therefore, the characteristic equation isand λ  = Re(λ) ± i Im(λ), where

Obviously, (i) Re(λ)∣(ϵ 2  = 0) = 0, ; (ii) ; (iii) from Lemma 3.3, we have W 1  < 0.

Therefore, by the Hopf bifurcation theory, we obtain

There exists a σ  >  0, and a function ϵ   =   ϵ (x ) defined on , which satisfies . And when ϵ   =   ϵ (x )  >  0, system (3.7) has a unique stable limit cycle which passes through (x ,   47/19).

The bifurcation diagram of system (3.5), (3.7) is in Fig. 3.1 . The bifurcation diagram shows that the unique positive equilibrium (x ∗,  y ∗) is unstable when 3738/361 <  A  < 12 and is stable when A  < 3738/361 or A  > 12. Thus Hopf bifurcations occur and there exists a stable limit cycle when 3738/361 <  A  < 12, which are the conclusions in Theorem 3.2, Theorem 3.4. As an example, setting A  = 11, we obtain the phase portrait of system (3.5) or (3.7) (see Fig. 3.2 ).

Fig. 3.1

The bifurcation diagram of (3.5) (or (3.7)) by parameter A and the y-coordinate of the limit cycle.

Fig. 3.2

A stable limit cycle of (3.5) when m = −3, n = 10, A = 11, q = 47, p = 9.

In the case that Δ > 0, R0 < 1 and 1 + q + pm < 0

From Theorem 1.3, we also know that when Δ  > 0, R 0  < 1 and 1 +  q  +  pm  < 0, system (1.3) has two positive equilibria E 1(x 1,  y 1) and E 2(x 2,  y 2), and E 1(x 1,  y 1) is a saddle. Furthermore, E 2(x 2,  y 2) is a center when S 2(x 2) = 0. The set of numbersatisfy thatNow we discuss the local property of the unique positive equilibrium E 2(x 2,  y 2) when (m,  n,  p,  q,  A) = (−5, 18, 7, 11, 6). In this case, we have x 2  = 1/9, y 2  = 11/9. We considerBy the similar argument as above, we obtain that the first Liapounov constant of system (3.8) about (1/9, 11/9) is . Therefore, we get

(1/9,   11/9) is a unstable weak focus of multiple one of system (3.8) .

In the following, we choose A as a bifurcation parameter. LetThen we may write (3.8) as follows:The second positive equilibrium of system (3.9) is and the linearizing matrix of (3.9) at E 2 iswhere

Therefore, the characteristic equation isand λ  = Re(λ) ± i Im(λ), where

Obviously, (i) Re(λ)∣(δ  = 0) = 0, ; (ii) ; (iii) from Lemma 3.5, we have W 1  > 0. Thus, by the Hopf bifurcation theory, we obtain the following conclusion.

There exists a σ  >  0, and a function defined on , which satisfies . And when system (3.9) has a unique unstable limit cycle which passes through .

Fig. 3.3 show the bifurcation diagram of system (3.9). When δ  ⩽ 0, the positive equilibrium E 2(x 2,  y 2) is unstable focus; when δ  > 0, the stability of the equilibrium E 2(x 2,  y 2) is changed and a Hopf bifurcation occur. Hence, when δ  > 0 the positive equilibrium E 2(x 2,  y 2) is a stable focus and there exists a unstable limit cycle. As an example, setting δ  = 0.008702, we obtain the phase portrait of system (3.9) (see Fig. 3.4 ).

Fig. 3.3

The bifurcation diagram of (3.9) by parameter δ and the y-coordinate of the limit cycle.

Fig. 3.4

A unstable limit cycle when m = −5, n = 18, A = 6.008702, q = 11, p = 7.

Conclusion

Epidemic mathematical models have become important tools to study the transmission dynamics of infectious diseases in host populations. There have been lots of works on the stability of endemic equilibrium for some epidemic models (cf. [10], [11], [12] and references therein). In this paper, we have studied the bifurcation behaviors of the SIR epidemic model with the general non-monotone and nonlinear incidence rate kIS/(1 +  βI  +  αI 2). From the analysis, we have found that there exist some values of the model such that the model can undergo a series of bifurcations, such as saddle-node bifurcation, Hopf bifurcation and Bogdanov–Takens bifurcation. When a stable limit cycle surrounds the endemic equilibrium, it means that the number of the infective tends to a periodic function and the disease will exhibit frequently regular oscillation. Hence, the disease become periodic outbreak as time evolves. On the other hand, there exist some parameters values such that the model has two endemic equilibria (one is a saddle and the other is center-type equilibrium) and a stable homoclinic loop. Thus, the disease will persist.