Complex dynamics in an SIS epidemic model induced by nonlinear incidence

We study an epidemic model with nonlinear incidence rate, describing the saturated mass action as well as the psychological effect of certain serious diseases on the community. Firstly, the existence and local stability of disease-free and endemic equilibria are investigated. Then, we prove the occurrence of backward bifurcation, saddle-node bifurcation, Hopf bifurcation and cusp type Bogdanov-Takens bifurcation of codimension 3. Finally, numerical simulations, including one limit cycle, two limit cycles, unstable homoclinic loop and many other phase portraits are presented. These results show that the psychological effect of diseases and the behavior change of the susceptible individuals may affect the final spread level of an epidemic.


Introduction
In generic, the population, for an SIS model, is always separated into two compartments, susceptible and infective individuals. And, in most SIS epidemic models (see Anderson and May [2]), the incidence takes the mass-action form with bilinear interactions. However, in practical, to describe the transmission process more realistic, it is necessary to introduce the nonlinear contact rates [1].
Actually, various forms of nonlinear incidence rates have been proposed recently. For example, in order to incorporate the effect of behavioral changes, Liu, Levin, and Iwasa [3] used a nonlinear incidence rate of the form where κI ι represents the force of infection of the disease, 1/(1 + αI h ) is a description of the suppression effect from the behavioral change of susceptible individuals when the infective population increases. ι, h and κ are all positive constants, and α is a nonnegative constant. See also Hethcote and van den Driessche [4], Moghadas [5] and Alexander and Moghadas [6,7], etc.
see Fig. 1(a), which can describe the effects of psychology function caused by protection measures and intervention policies. Obviously, g 1 (I) is increasing with small I and decreasing with large I, that is, g 1 (I) is non-monotone. The function g 1 (I) can be used to interpret the psychological effect: for a very large number of infective individuals, the infection force may decrease as the number of infective individuals increases, since large number of infective may lead to the reducing of the number of contacts per unit time. For example, in 2003, epidemic outbreak of severe acute respiratory syndrome (SARS) had such psychological effects on the general public (see [11]), and aggressive measures and policies had been taken, such as border screening, mask wearing, quarantine, isolation, etc. And also [11] showed that either the number of infective individuals tends to zero as time evolves or the disease persists. Furthermore, Li, Zhao and Zhu (see [20]) studied the following SIS model, which describes behavior change effect of susceptible individual when infectious population increases dS dt = Λ − d 0 S − g 2 (I)S + δI, dI dt = g 2 (I)S − (d 0 + d + δ)I, where g 2 (I) = αI 2 1 + β 0 I 2 , see Fig. 1(b). By the qualitative and bifurcation analyses, they showed that the maximal multiplicity of weak focus is 2, and proved that the model can undergo a Bogdanov-Takens bifurcation of codimension 2. These results illustrate that the behavior change of the susceptible individuals may affect the final spread level of an epidemic.
Actually, both the effect of psychological and behavior affect of susceptible individuals have influence on the transmission of the disease. Thus, motivated by the above researches, we consider a nonlinear incidence rate of a SIS model as follows where G(I) = aI 2 c + I 2 + bI c + I 2 , see Fig. 1(c), which can describe the effect of psychology and behavior change of susceptible individuals. S(t) and I(t) represent the number of susceptible individuals and infected individuals at time t, respectively. Λ is the recruitment rate of population, d is the natural death, µ is the disease-induced death rate, and σ represents the recovered rate, a, b and c are all positive constants.
The organization of this paper is as follows. First of all, we analyze the existence of equilibria and local stability of equilibria. Then, we study the existence of Hopf bifurcation around the positive equilibrium at the critical value under the conditions of R 0 < 1 and R 0 > 1. We also show that these positive equilibria can be weak focus for some parameter values and a cusp type of Bogdanov-Takens bifurcation of codimension 3. Finally, we give some brief conclusion.

Types and Stability of equilibria
In this section, our aim is to perform an elaborative analysis of equilibria for system (2).
Proof Summing up the two equations in (3), we can get Thus, lim sup t→∞ (x + y) ≤ Λ d , which implies the conclusion. Obviously, system (3) always has a unique disease-free equilibrium E 0 = ( Λ d , 0). The positive equilibria of (3) can be gained by solving the following algebraic equations Λ − xy(ay + 1) y 2 + c − dx + σy = 0, xy(ay + 1) Denote the basic reproduction number as follows, Computing the discriminant of (4), we get Furthermore, define then ∆ > 0 if and only if R 0 > R * , ∆ = 0 if and only if R 0 = R * and ∆ < 0 if and only if R 0 < R * . It is obvious that R * < 1 and we can obtain the following theorem.
Theorem 1 Model (3) always has a disease-free equilibrium E 0 and the following conclusions hold.
(i) When R 0 < 1, we have (a) if R 0 < R * , then system (3) has no positive equilibrium; (b) if R 0 = R * and Λ > (1 − σ)/a, then system (3) has a unique positive equilibrium ay k +1 (k = 2, 3) and In the following, we discuss the local stability of E k (x k , y k ) (k = 0, 1, 2, 3, 4, 5) and present the corresponding phase portrait. By directly calculating, the Jacobian matrix at equilibrium E k is We have In addition, after some complicated computations, we get .
Indeed, if R 0 = 1, the Jacobian J 0 is diagonalizable with eigenvalues λ 1 = −d and λ 2 = 0 and respective eigenvectors ν 1 = (1, 0) and ν 2 = (− 1−σ d , 1). By the transformation u system (5) can be rewritten as If 2a − 1−σ dc = 0, according to the calculation of center manifold [15], we know that the center manifold x = h(y) of (5) begins with quadratic term of y. In addition, from the second equation of (5), we can easily obtain that the equation restricted to the center manifold as follows By applying the Theorem 7.1 in Zhang et al. [15], E 0 is a saddle-node. If 2a − 1−σ dc = 0, then the center manifold turns into dy dt = 1 c y 3 + O(y 4 ).
Since 1 c > 0 and the first nonzero item is uneven. Thus, the equilibrium E 0 is a repelling node, according to Theorem 7.1 in Zhang et al. [15].
From the expression of ψ 1 and φ 1 , we can see that one of the eigenvalues of characteristic matrix of E 1 is zero and the other is nonzero if ψ 1 = 0. The type of E 1 can be directed proved by checking the conditions in Zhang et al. ( [15], Theorem 7.1-7.3). So, we have the following results.
Theorem 4 Suppose that R * < R 0 < 1 and aΛ > 1 − σ, then system (3) has two positive equilibria E 2 and E 3 , and equilibrium E 2 is a hyperbolic saddle for all permissible choices of parameters, equilibrium E 3 is not degenerate. Moreover, Proof Note that φ 2 is less than zero since ∆ < (aΛ + 1 − σ) 2 , then E 2 is a hyperbolic saddle for any choices of parameters. And at E 3 , we have φ 3 > 0. Thus, the stability of equilibrium E 3 depends on the sign of ψ 3 .
When ψ i = 0 (i = 3, 5), the dynamics of system (3) can be easily seen in Fig. 2, Fig.  3 and Fig. 4, respectively. The dynamical behaviors of system (3) when ψ i = 0 (i = 3, 5) will be discussed in detail in the next section. Remark 1 In fact, Fig. 2 (a) shows the occurrence of bi-stability, in which solution may converge to one of the two equilibria, depending on the initial conditions. And in practical, this interesting phenomenon implies that initial states determine whether the disease can dies out or not.
Remark 2 From Fig. 2 (a), we can see there exists two seperatrices. All solutions tend to the disease-free equilibrium E 0 except the two green lines tend to equilibrium E 2 .
Theorem 7 Suppose that R 0 = 1 and Λ > 1−σ a , then system (3) has a unique positive equilibrium E 4 . If ψ 4 > 0, then there exists at least one stable limit cycle in the interior of the first quadrant.
Proof Indeed, the Jacobian J 0 is diagonalizable with eigenvalues λ 1 = −d and λ 2 = 0 and respective eigenvectors ν 1 = (1, 0) and ν 2 = (− 1−σ d , 1). By the transformation  system (5) becomes The theorem of Chochitaichvili [21] yields directly that system (6) is topologically equivalent to the system: It is easy to see that 2a − 1−σ dc > 0 according to the condition of R 0 = 1 and Λ > 1−σ a , and then we can gain that there exist a unique repelling equilibrium E 4 in the region D 1 shown in (a) of Fig.5. Consequently, by the Poincaré-Bendixson theorem, at least one stable limit cycle appears in the interior of the first quadrant. Similarly, when ψ 5 > 0, we have the following result.
Theorem 8 Suppose that R 0 > 1. If ψ 5 > 0, then there exists at least one stable limit cycle in the interior of the first quadrant.
Proof Indeed, the Jacobian J 5 has eigenvalues λ 1 = −d and λ 2 = R 0 − 1 > 0, when R 0 > 1. Thus, we can gain that there exist a E 5 , which is the unique repelling equilibrium in the region D 2 shown in (b) of Fig.5. Consequently, by the Poincaré-Bendixson theorem, at least one stable limit cycle appears in the interior of the first quadrant.
a and ψ 4 > 0, there exists a stable limit cycle. (b): When R 0 > 1 and ψ 5 > 0, there exists a stable limit cycle. Fig. 3, we can see clearly that there exist a stable limit cycle enclosing the equilibrium E 5 even though E 0 is a saddle-node. Similarly, from (b) of Fig. 4, we can also gain that there exists a stable limit cycle if E 4 is unstable.

Bifurcations Backward bifurcation
Remark 4 When R 0 = 1 and Λ > 1−σ a , system (3) exhibits a unique positive equilibrium E 4 , which means that once R 0 crosses 1, the disease can invade to a relatively high level. And this is one of the main characters of backward bifurcation [19].
Remark 5 Actually, backward bifurcation did not emerge with a = 0, which is considered in [20]. This indicates that introducing the non-motonic incidence into model (1) makes the epidemic model more complex and exhibit richer dynamical behaviors.

Hopf bifurcation
In this subsection, we will study the Hopf bifurcation of system (3) for (i) R * < R 0 < 1 and Λ > (1 − σ)/a; (ii) R 0 > 1. From the discussion in Section 2, it can be seen that Hopf bifurcation may occur at E 3 , E 5 . As the expressions of equilibria E 3 and E 5 are same, no considering the values of every parameters. Based on Theorem 4, Theorem 5 and Theorem 6, we know that the stability of E 3 and E 5 are similar and when ψ k = 0, E k (k = 3, 5) is a weak focus or center. Thus, we show the existence of Hopf bifurcation around E k (k = 3, 5).
Theorem 10 Suppose E k (k = 3, 5) exist, then model (2) undergoes a Hopf bifurcation at equilibrium E k if ψ k = 0. Moreover, (a) if η < 0, there is a family of stable periodic orbits of model (3) as ψ k decreases from 0; (b) if η = 0, there are at least two limit cycles in (3), where η will be defined below; (c) if η > 0, there is a family of unstable periodic orbits of (3) as ψ k increases from 0.
Proof From the above discussions, we can see that trJ k = 0 if and only if ψ k = 0, and det J k > 0 when equilibrium E k exists. Therefore, the eigenvalues of J k are a pair of pure imaginary roots if ψ k = 0. From direct calculations we have that By Theorem 3.4.2 in [16], Ψ k = 0 is the Hopf bifurcation point for (3). Next, introduce a new time variable τ by dt = (y 2 + c)dτ . By rewriting τ as t, we obtain the following equivalent system of (3) dy dt = xy(ay + 1) − y(y 2 + c).
Let X = x − x k and Y = y − y k , still use (x, y) to express (X, Y ), then system (7) becomes LetÊ denote the origin of x − y plane. Since E k satisfies Eq. (3), we obtain and it is easy to verify that b 11 + b 22 = 0 if and only if ψ k = 0. Let ω = (det J(Ê)) 1 2 , u = −x and v = b11 ω x + b12 ω y, then the normal form of system (7) reads where f uv denotes ∂ 2 f ∂u∂v (0, 0), etc. Then, by computing we obtain By Theorem 3.4.2 and Theorem 3.4.11 in [16], the rest claims in Theorem 10 are true.
Remark 6 What need to be noted here is that the expression b 12 in Theorem 10 is nonzero. Otherwise, we have det J k = b 11 b 22 < 0, since b 11 + b 22 = 0, which is an contradiction.
Next, we present examples to show that equilibrium E k can be a stable weak focus of multiplicity two, and under a small perturbation, system (3) undergoes a degenerate Hopf bifurcation and produces two limit cycles.
With the aid of Mathematica, we get Besides, we give the numerical simulation graphs for one limit cycle and two limit cycles under small perturbations of some parameters. From Fig. 6 (a) , we can see that there exist only one limit cycle around the endemic E 3 , and Fig. 6 (b) shows us that a new limit cycle emerge with small perturbations of parameters Λ and c. It is worth emphasizing that if we change the values of parameters Λ and c, an unstable homoclinic loop arise, which is shown in Fig.7.
Around equilibrium E 5 , we can gain the same result from Fig. 8. Under the condition that parameter a, d, σ take value 5.6923, 0.11, 0.2 and change value of Λ and c from 0.45102, 3.6062 to 0.51, 3.5972, respectively, which is a minor change, then the number of limit cycles will add one. Thus, the interior equilibrium E 5 is a stable weak focus of multiplicity two if (σ, Λ, d, a, c) = (0.2, 0.51, 0.11, 5.6923, 3.5962).

Remark 7
As a matter of fact, the reproduction number is equal to zero in [8,20], which simplifies the condition that Hopf bifurcation occur. In our model, we also comprehensively discuss the existence of Hopf bifurcation when R 0 < 1, R 0 = 1 and R 0 > 1. Besides, authors in [20] did not show the appearance of homoclinic loop, which is an interesting bifurcation phenomenon given in Fig.7.
is equivalent to the system in some small neighborhood of (0, 0) after changes of coordinates.

Lemma 4 The system
is equivalent to the system in some small neighborhood of (0, 0) after changes of coordinates.
From the proving process of Theorem 4, it is easily to know that detJ 1 = 0. And Theorem 3 suggests that the characteristic matrix J 1 possesses a zero eigenvalue with multiplicity 2 when ψ 1 = 0, which shows that system (3) may admit a Bogdanov-Takens bifurcation. Thus, we can prove the following theorem.
Define two functions: Theorem 11 Suppose that R 0 = R * and ψ 1 = 0, then the only interior equilibrium E 1 of system (3) is a cusp. Moreover, (a) if f (y 1 )g(y 1 ) = 0, then E 1 is a cusp of codimension 2; (b) if f (y 1 )g(y 1 ) = 0, then E 1 is a cusp of codimension greater than or equal to 3.
Proof Changing the variables as X = x − x 1 , Y = y − y 1 , system (3) becomes where Make the non-singular linear transformation 1 0 x y .
Unfortunately, because of the complexity of f (y 1 ) and g(y 1 ), we can not determine which of these three situations can occur theoretically. But we will show for some parameter values, f (y 1 ) = 0 and g(y 1 ) = 0, i.e. E 1 is a cusp point. In the following, we will give an example to show that the case (b) of Theorem 11 occur.

Remark 10
Xiao and Ruan (see [10]) showed that either the number of infective individuals tends to zero as time evolves or the disease persists. Authors in [8,20] proved that their epidemic model with saturated incidence rate, undergo a Bogdanov-Takens bifurcation of codimension 2. When we consider the incidence of a combination of the saturated incidence rate and a non-monotonic incidence, the codimension of Bogdanov-Takens bifurcation can be up to 3.

Conclusions
In this paper, by combing qualitative and bifurcation analyses, we study an SIS epidemic model with the incidence rate aI 2 c+I 2 + bI c+I 2 , which can describe the inhibition effect from the behavioral change and interpret the "psychological" effect. In Section 2, we give a full-scale analysis for the types and stability of equilibria E i (i = 0, 1, 2, 3, 4, 5). We prove that system (2) can occur backward bifurcation and the backward bifurcation will disappear if a = 0. At equilibrium E i (i = 3, 5), degenerate Hopf bifurcation arises under certain conditions. When the critical condition Ψ i (i = 3, 5) satisfied, we calculate the Liapunov value of the weak focus and obtain the value 2 for the maximal multiplicity of weak focus, indicating that there exist at most two limit cycles around E i (i = 3, 5). In Fig. 6, and Fig. 8, we give the phase portrait corresponding to equilibrium E 3 and E 5 about exhibiting a unique limit cycle and adding a new limit cycle after small perturbation of parameter Λ and c. In Subsection 3.3, we proved that the model exhibits Bogdanov-Takens bifurcation of codimension 2 and codimension 3, under certain conditions. If parameter a = 0, the model can just have Bogdanov-Takens bifurcation of codimension 2, shown in [20].
In fact, we show that the model exhibits multi-stable states. This interesting phenomenon indicates that the initial states of an epidemic can determine the final states of an epidemic to extinct or not. Moreover, the periodical oscillation signify that the trend of the disease may be affected by the behavior of susceptible and the effect of psychology of the disease.