key: cord-1009049-7vs4zk3c
authors: Lu, Min; Huang, Jicai; Ruan, Shigui; Yu, Pei
title: Bifurcation analysis of an SIRS epidemic model with a generalized nonmonotone and saturated incidence rate()
date: 2019-07-15
journal: J Differ Equ
DOI: 10.1016/j.jde.2019.03.005
sha: 96bbfb4e9db97a3c63b247ef0898122f10d099dc
doc_id: 1009049
cord_uid: 7vs4zk3c

In this paper, we study a susceptible-infectious-recovered (SIRS) epidemic model with a generalized nonmonotone and saturated incidence rate [Formula: see text] , in which the infection function first increases to a maximum when a new infectious disease emerges, then decreases due to psychological effect, and eventually tends to a saturation level due to crowding effect. It is shown that there are a weak focus of multiplicity at most two and a cusp of codimension at most two for various parameter values, and the model undergoes saddle-node bifurcation, Bogdanov-Takens bifurcation of codimension two, Hopf bifurcation, and degenerate Hopf bifurcation of codimension two as the parameters vary. It is shown that there exists a critical value [Formula: see text] for the psychological effect, and two critical values [Formula: see text] for the infection rate such that: (i) when [Formula: see text] , or [Formula: see text] and [Formula: see text] , the disease will die out for all positive initial populations; (ii) when [Formula: see text] and [Formula: see text] , the disease will die out for almost all positive initial populations; (iii) when [Formula: see text] and [Formula: see text] , the disease will persist in the form of a positive coexistent steady state for some positive initial populations; and (iv) when [Formula: see text] and [Formula: see text] , the disease will persist in the form of multiple positive periodic coexistent oscillations and coexistent steady states for some positive initial populations. Numerical simulations, including the existence of one or two limit cycles and data-fitting of the influenza data in Mainland China, are presented to illustrate the theoretical results.

In the classical infectious disease transmission models (Kermack and McKendrick [16] , Hethcote [12] ), the population is divided into three classes labeled by S(t), I (t) and R(t), representing the numbers of susceptible, infective, and recovered or removed individuals at time t , respectively. Assuming the recovered individuals have temporary immunity, the classical susceptibleinfectious-recovered (SIRS) model can be written as follows:

where b is the recruitment rate of the population, d is the natural death rate of the population, μ is the natural recovery rate of the infective individuals, δ is the rate at which recovered individuals lose immunity and return to the susceptible class. g(I )S is called the incidence rate, and g(I ) is a function to measure the infection force of a disease. The incidence rate g(I )S plays a very important role in describing the evolution of infectious disease. In [16] Kermack and McKendrick assumed that the incidence rate is bilinear, i.e., g(I )S = kI S, (1.2) where g(I ) = kI is unbounded when I ≥ 0. The bilinear incidence rate (1.2) might be true when the number of the infective individuals I (t) is small, but it becomes unrealistic when I (t) is getting larger. Studying the data on the cholera epidemic spread in Bari, Italy, in 1973, Capasso et al. [5] and Capasso and Serio [6] proposed two types of nonlinear incidence rates, namely a saturated incidence rate and an incidence rate taking psychological effect into account (see Fig. 1.1) .

(a) Saturated Incidence Rates. To compare with the bilinear incidence rate (1.2), Capasso and Serio [6] proposed a saturated incidence rate of the following form:

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 or from the crowding effect of the infective individuals (Ruan and Wang [25] ). Notice that the nonlinear incidence function g(I ) eventually tends to a saturation level k α when I is getting larger. The general incidence rate g(I )S = kI p S 1 + αI q (1.4) was proposed by Liu et al. [20] (q = p − 1) and Hethcote and van den Driessche [13] (p ≥ q), and used by a number of authors, see, for example, Alexander and Moghadas [1, 2] , Derrick and van den Driessche [9, 10] , Li et al. [17] , Liu et al. [19] , Lizana and Rivero [21] , Moghadas and Alexande [22] , and Wang [28] . According to Tang et al. [27] and Hu et al. [14] , the incidence function g(I ) = kI p 1+αI q includes three types: (i) unbounded incidence function: p > q; (ii) saturated incidence function: p = q; and (iii) nonmonotone incidence function: p < q. In order to better understand the generic bifurcations in SIRS models with saturated incidence rates, Ruan and Wang [25] studied model (1.1) with a specific nonlinear incidence rate, g(I )S = kI 2 S 1 + αI 2 , (1.5) and presented a detailed qualitative and bifurcation analysis of the model. (b) Nonmonotone Incidence Rates with Psychological Effect. Though Capasso and Serio [6] described incidence functions taking into account psychological effect, but they did not give any specific function to describe such incidence. To model the effects of psychological factor, protection measures and intervention policies when a serious disease emerges, Xiao and Ruan [29] proposed the following specific incidence rate:

where the incidence function g(I ) = kI 1+αI 2 is nonmonotone when I ≥ 0 (see Fig. 1.1(b) ). This implies that the contact rate and the infection probability are increasing when a new infectious disease emerges, since people have very little knowledge about the disease. However, when I is large and the disease becomes more serious, psychological factor leads people to implementing measures to control the spread of the disease. For example, in the outbreak of severe acute respiratory syndrome (SARS), the aggressive measures and policies, such as border screening, mask wearing, quarantine, isolation, etc. were proved to be very effective in reducing the transmission. So the infection force decreases when the number of infected individuals becomes larger. Liu et al. [18] used (1.6) to describe the psychological effect toward avian influenza in the human population. In [29] , Xiao and Ruan presented a global analysis for model (1.1) with nonmonotone incidence rate (1.6) , and showed that either the number of infective individuals tends to zero as time evolves or the disease persists. Hence, model (1.1) with nonmonotone incidence rate (1.6) cannot exhibit complicated dynamics and bifurcations.

Xiao and Zhou [30] considered a complete form of the nonmonotone incidence rate (1.6) as follows:

where β is a parameter such that 1 + βI + αI 2 > 0 for all I ≥ 0, hence, β > −2 √ α. They presented qualitative analysis of model (1.1) with nonmonotone incidence rate (1.7) and showed the existence of bistable phenomenon and periodic oscillation. Zhou et al. [32] further studied the existence of different kinds of bifurcations, such as Hopf and Bogdanov-Takens bifurcations.

(c) Nonmonotone and Saturated Incidence Rates. Note that both nonmonotone incidence functions kI 1+αI 2 and kI 1+βI +αI 2 , given in (1.6) and (1.7), tend to zero when I goes to infinite, which indicates that the psychological or inhibitory effect from the behavioral change of the susceptible individuals or from the crowding effect of the infective individuals is too strong, which may be unreasonable for some specific infectious diseases, such as influenza. Thus, a more reasonable incidence function g(I ) may be one that first increases to a maximum when a new infectious disease emerges or an old infectious disease reemerges, then decreases due to psychological effect, and eventually tends to a saturation level due to crowding effect. On the other hand, in some specific infectious diseases, the incidence rate may not be monotonic or non-monotonic alone, a more general incidence rate may have a combination of monotonicity, nonmonotonicity and saturation properties.

Based on the above discussions, in this paper, we propose the following general nonmonotone and saturated incidence rate:

where α > 0 is a parameter which measures the psychological or inhibitory effect, and k is the infection rate. β is a parameter such that 1 + βI + αI 2 > 0 for all I ≥ 0, hence, β > −2 √ α. When β ≥ 0, g(I ) = kI 2 1+βI +αI 2 is monotonic, always increases and then tends to a saturated level k α as I goes to infinite. When −2 √ α < β < 0, g(I ) is nonmonotonic, increases when I is small and decreases when I is large, and finally tends to a saturated level k α as I goes to infinite (see Fig. 1 .2). It can be seen from Fig. 1 .2 that when α and k are fixed, g(I ) increases faster for smaller β when I is small, and then decreases to the same fixed value k α for any β. When β = 0, the general incidence rate (1.8) becomes the saturated incidence rate (1.5). Ruan and Wang [25] studied model (1.1) with incidence rate (1.5) and presented a detailed qualitative and bifurcation analysis of the model. They showed that the system undergoes Bogdanov-Takens bifurcation of codimension two and Hopf bifurcation of codimension one, but the exact codimension of Hopf bifurcation remains unknown. In [27] , Tang et al. proved that the maximal multiplicity of the weak focus is two by complex algebraic calculations but did not give the explicit conditions for degenerate Hopf bifurcation. Moreover, they claimed the existence of Bogdanov-Takens bifurcation of codimension three. They also conjectured that model (1.1) with nonmonotone incidence rate would not exhibit complicated dynamics and bifurcations.

In this paper, for the case β > −2 √ α we find that model (1.1) with general incidence rate (1.8) can exhibit complicated dynamical behaviors and bifurcation phenomena. More precisely, we will show that there are a weak focus of multiplicity at most two and a cusp of codimension at most two for various parameter values, and the model undergoes saddle-node bifurcation, Bogdanov-Takens bifurcation of codimension two, Hopf bifurcation, and degenerate Hopf bifurcation of codimension two as the parameters vary. By considering a more general model, our results show that Bogdanov-Takens bifurcation of codimension higher than two cannot occur, which clarifies and corrects the corresponding results obtained by Tang et al. [27] , who claimed the existence of Bogdanov-Takens bifurcation of codimension three. Moreover, our results about degenerate Hopf bifurcation of codimension two can be seen as a complement of the results obtained by Ruan and Wang [25] , who only discussed the Hopf bifurcation of codimension one. Furthermore, it is shown that there exists a critical value α = α 0 for the psychological or inhibitory effect, and two critical values k = k 0 , k 1 (k 0 < k 1 ) for the infection rate such that: (i) when α > α 0 , or α ≤ α 0 and k ≤ k 0 , the disease will die out for all positive initial populations; (ii) when α = α 0 and k 0 < k ≤ k 1 , the disease will die out for almost all positive initial populations; (iii) when α = α 0 and k > k 1 , the disease will persist in the form of a positive coexistent steady state for some positive initial populations; and (iv) when α < α 0 and k > k 0 , the disease will persist in the form of multiple positive periodic coexistent oscillations and coexistent steady states for some positive initial populations.

The rest of the paper is organized as follows. In section 2, we carry out a qualitative analysis to give the types and stability of equilibria of model (1.1) with general incidence rate (1.8) . In section 3, we show that the model undergoes saddle-node bifurcation, Bogdanov-Takens bifurcation of codimension at most two, Hopf bifurcation, and degenerate Hopf bifurcation of codimension at most two as parameters vary. In section 4 we use the model to simulate the influenza data in Mainland China from 2004 to 2017. Conclusion and discussion are given in section 5.

We consider an SIRS epidemic model in the following form

and R(t) denote the numbers of susceptible, infective, and recovered individuals at time t , respectively, b > 0 is the recruitment rate of the population, d > 0 is the natural death rate of the population, k > 0 is the infection rate, μ > 0 is the natural recovery rate of the infective individuals, δ > 0 is the rate at which recovered individuals lose immunity and return to the susceptible class, α > 0 is a parameter which measures the psychological or inhibitory effect, and β > −2 √ α such that 1 + βI + αI 2 > 0 for all I ≥ 0.

To study the dynamics of model (2.1), we first present a lemma.

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

Proof. Summing up the three equations in (2.1) and denoting N(t) = S(t) + I (t) + R(t), we have

It is clear that N(t) = b d is a solution and for any N(t 0 ) ≥ 0, the general solution is

Thus,

which implies the conclusion. 2 It is clear that the limit set of system (2.1) is on the plane S + I + R = b d . Thus, we focus on the reduced system

We know that the positively invariant set of system (2.2) is

System (2.2) always has an equilibrium E 0 = (0, 0) which corresponds to the disease-free equilibrium ( b d , 0, 0) of system (2.1). To find the endemic equilibria, we rescale (2.2) by using

Then we obtain (for simplicity we still denote τ by t)

where

It can be seen that (2.4) and the positively invariant region of system (2.3) is

is an invariant line, and when y = 0, we have dy dt > 0. To find the positive equilibria of system (2.3), we set

which yield

We can see from (2.5) and (2.6) that system (2. 3) has at most two positive equilibria E 1 (x 1 , y 1 ) and E 2 (x 2 , y 2 ), which may coalesce into a unique positive equilibrium E * (x * , y * ), where

The discriminant of (2.6) is

and we have

Then we have the following existence conditions of equilibria in system (2.3). where 0 < x * ≤ 1. (III) System (2.3) has two positive equilibria E 1 (x 1 , y 1 ) and E 2 (x 2 , y 2 ) if and only if n < n * and A > mp, where 0 < x 1 < x * < x 2 < 2x * .

Next we study local stability of equilibria of system (2.3). We first study the disease-free equilibrium E 0 (0, 0). The Jacobian matrix of system (2.3) at E 0 (0, 0) is 

which has two eigenvalues λ 1 = −1, λ 2 = −p. We obtain the following result. To discuss whether the disease can invade the population, we study the global stability of the equilibrium (0, 0). Since x = 0 is an invariant line and D is positively invariant, by index theory, we can conclude that system (2.3) does not have nontrivial periodic orbits in R 2 + when system (2.3) has no positive equilibria. 

while n ≤ n * and A ≤ mp are equivalent to α ≤ α 0 and k ≤ k 0 , where

These indicate that the disease will die out for all positive initial populations when the psychological or inhibitory effect α is greater than a critical value α 0 , or the psychological effect α is smaller than the critical value α 0 and the infection rate k is also smaller than a critical value k 0 . Now we consider the positive equilibria of system (2.3). The positive equilibria with coordinates (x, y) satisfy

and the Jacobian matrix of system (2.3) at a positive equilibrium E(x, y) is given by

and its sign is determined by

The trace of J (E) is

and its sign is determined by

To discuss the topological type of the positive equilibria of system (2.3), we let

which will be used in the following theorem. (I) if q = q 1 , then E * (x * , y * ) is a saddle-node, which is attracting (or repelling) if q < q 1 (or q > q 1 ); (II) if q = q 1 , then E * (x * , y * ) is a cusp of codimension two.

The phase portraits are shown in Fig. 2 

Proof. Substituting x * and n * into S D and S T , we deduce that S D (x * ) = 0 and

Then, using Theorem 7.1 in Zhang et al. [33] , we obtain the conclusion in (I).

To show that the assertion (II) holds, let X = x − x * , Y = y − y * , n = n * and q = q 1 . Then using Taylor expansions, we can rewrite system (2.3) as follows (for simplicity, we still denote X, Y by x, y, respectively)

where

(2.14)

By Remark 1 of section 2.13 in [24] (see also [15] ), we obtain an equivalent system of (2.14) in the small neighborhood of (0, 0) as follows:

It is seen that

and substituting it into n = n * and q = q 1 , we obtain

Further, since n > 0, p > q > 0 and m > −2 √ n, it follows from (2.16) that

. By the results in [15] , E * (x * , y * ) is a cusp of codimension two. 2 Remark 2.2. When m = 0 (i.e., β = 0), our system (2.2) is reduced to system (1.3) in [27] , our Theorem 2.5 (II) indicates that Bogdanov-Takens bifurcation of codimension higher than two cannot occur around E * , which clarifies and corrects the corresponding results in [27] , where the authors claimed the existence of Bogdanov-Takens bifurcation of codimension three.

When the psychological effect α equals the critical value α 0 given by (2.8), Theorem 2.5 implies that whether the disease persists or dies out will depend on the infection rate and the initial population. More precisely, when the psychological effect α equals the critical value α 0 and the infection rate k is greater than the first critical value k 0 given by (2.9), i.e., n = n * and A > mp, then system (2.1) has two equilibria, a disease-free equilibrium and an endemic equilibrium. Moreover, from Fig. 2 .2 (b) and (c), we can see that the disease will die out for almost all positive initial populations if the infection rate k is less than another larger critical value (2.17) and the disease will persist to a positive coexistent steady state for some positive initial population if the infection rate k is greater than the second critical value k 1 (k > k 1 , i.e. q < q 1 ), as shown in Fig. 2 

Theorem 2.6. When n < n * , A > mp and conditions (2.4) are satisfied, then system (2.3) has two positive equilibria E 1 (x 1 , y 1 ) and E 2 (x 2 , y 2 ). Moreover, E 1 is always a hyperbolic saddle and E 2 is

where S T is given in (2.11).

To determine the types of E 1 and E 2 , it is suffice to consider the signs of S D (x 1 ), S D (x 2 ) and S T (x 2 ), where S D and S T are given in (2.10) and (2.11), respectively, from which we have

On the other hand, since x 1 and x 2 are two different positive roots of (2.6), we have

yielding S D (x 1 ) < 0 and S D (x 2 ) > 0. Hence we obtain the types of E 1 and E 2 . 2 Remark 2.4. It is seen from Theorem 2.6 that when the psychological effect α is smaller than the critical value α 0 but the infection rate k is greater than the first critical value k 0 , i.e., n < n * and A > mp, then the disease will become more severe because of the existence of multiple positive coexistent steady states.

From Theorems 2.2, 2.5 and 2.6, we know that system (2.3) may exhibit saddle-node bifurcation, Bogdanov-Takens bifurcation around the equilibrium E * (x * , y * ), and Hopf bifurcation around the equilibrium E 2 (x 2 , y 2 ). In this section, we investigate various possible bifurcations in system (2.3).

From Theorems 2.2, 2.5 and 2.6, we know that the surface

is the saddle-node bifurcation surface. When the parameters are varied to cross the surface from one side to the other side, the number of positive equilibria of system (2.3) changes from zero to two, the saddle-node bifurcation yields two positive equilibria. This implies that there exists a critical psychological effect value α 0 such that the disease cannot invade the population when α > α 0 (i.e., n > n * ), and the disease will persist for some positive initial populations when α ≤ α 0 (i.e., n ≤ n * ).

In this subsection, we discuss if system (2.3) undergoes Bogdanov-Takens Bifurcation of codimension two under a small parameter perturbation if the bifurcation parameters are chosen suitably. Actually, we have the following theorem. Bogdanov-Takens singularity 

where λ 1 and λ 2 are parameters in a small neighborhood of (0, 0). We are interested only in the phase portraits of system (3.1) when x and y lie in a small neighborhood of the interior equilibrium E * (x * , y * ). Let X = x − x * , Y = y − y * . Then system (3.1) can be rewritten as (for simplicity, we still denote X, Y by x, y, respectively)

where P 1 (x, y, λ 1 , λ 2 ) is a C ∞ function at least of third order with respect to (x, y), whose coefficients depend smoothly on λ 1 and λ 2 , and

and rewrite X, Y as x, y, respectively. Then system (3.2) becomes dx dt = y,

where Q 1 (x, y, λ 1 , λ 2 ) is a C ∞ function at least of third order with respect to (x, y), whose coefficients depend smoothly on λ 1 and λ 2 , and

, and rewriting X, Y as x, y, respectively, we obtain

where Q 2 (x, y, λ 1 , λ 2 ) is a C ∞ function at least of third order with respect to (x, y), whose coefficients depend smoothly on λ 1 and λ 2 , and

Note that when λ 1 = λ 2 = 0,

Further, let

Then system (3.5) can be rewritten as (still denote X, Y by x, y, respectively)

where Q 3 (x, y, λ 1 , λ 2 ) is a C ∞ function at least of third order with respect to (x, y), whose coefficients depend smoothly on λ 1 and λ 2 , and

Making the final change of variables by setting

then we finally obtain (still denote X, Y and τ by x, y and t, respectively)

where Q 4 (x, y, λ 1 , λ 2 ) is a C ∞ function at least of third order with respect to (x, y), whose coefficients depend smoothly on λ 1 and λ 2 , and

We can express μ 1 and μ 2 in terms of λ 1 and λ 2 as follows:

for m = A(p−1) p 2 , p > 0, q 1 > 0 and A > mp, the parameter transformation (3.8) is a homeomorphism in a small neighborhood of the origin, and μ 1 and μ 2 are independent parameters.

The results in Bogdanov [3, 4] and Takens [26] now imply that system (3.7) (i.e., (3.1) or (2.3)) undergoes Bogdanov-Takens bifurcation when (λ 1 , λ 2 ) changes in a small neighborhood of (0, 0). 2

By the results of Perko [24] , we obtain the following local representations of the bifurcation curves up to second-order approximations:

(i) The saddle-node bifurcation curve is

(ii) The Hopf bifurcation curve is

(iii) The homoclinic bifurcation curve is

The Bogdanov-Takens bifurcation diagram and corresponding phase portraits of system (3.1) are given in Fig. 3 .1, where we fix A = 2, m = −1, p = 2, q = 1/3 and n = 1/3, which satisfy n = n * , q = q 1 , A > mp and conditions (2.4), and then perturb A and q. The bifurcation curves H, HL and SN divide the small neighborhood of the origin in the parameter (λ 1 , λ 2 )-plane into four regions (see Fig. 3.1(a) ).

(a) When (λ 1 , λ 2 ) = (0, 0), the unique positive equilibrium is a cusp of codimension two (see Fig. 2.2(c) ).

(b) There are no equilibria when the parameters lie in region I (see Fig. 3.1(b) ), implying that the disease dies out.

(c) When the parameters lie on the curve SN, there is a unique positive equilibrium E * , which is a saddle-node. (d) Two positive equilibria, one is an unstable focus E 2 and the other is a saddle E 1 , will occur through the saddle-node bifurcation when the parameters cross the SN curve into region II (see Fig. 3.1(c) ).

(e) An unstable limit cycle will appear through the subcritical Hopf bifurcation around E 2 when the parameters are varied to cross the H curve into region III (see Fig. 3.1(d) ), where the focus E 2 is stable, whereas the focus E 2 is an unstable one with multiplicity one when the parameters lie on the curve H.

(f) An unstable homoclinic orbit will occur through the homoclinic bifurcation around E 1 when the parameters pass region III and lie on the curve HL (see Fig. 3.1(e) ).

(g) The relative location of one stable and one unstable manifold of the saddle E 1 will be reversed when the parameters are varied to cross the HL curve into region IV (compare Fig. 3.1(c) and Fig. 3.1(f) ).

In this subsection we discuss Hopf bifurcation around the equilibrium E 2 (x 2 , y 2 ). Let

which will be used in the following subsection.

To simplify the computation, by following the technique in [11] , we reduce system (2.3) to an equivalent polynomial differential system by using the following state variable scaling and time rescaling: 

(3.11)

Then taking the parameter scaling,

into system (3.11) and dropping the bars, we obtain

(3.12)

Since system (3.12) has an equilibrium Ẽ 2 (1, 1) (i.e., E 2 (x 2 , y 2 ) of system (2.3)), we have

Note that x 1 x 2 = p np+q+1 for system (2.3), which becomes x 1 x 2 = p np+q+1 under the parameter scaling (3.10), dropping the bars, so we have np + q + 1 − p > 0 because x 1 x 2 < 1. Regarding the above parameter scaling, the condition (2.4), n < n * and A > mp become

Since the transformation (3.10) is a linear sign-reserving transformation, system (3.12) and system (2.3) have the same qualitative property. Next letting dt = (1 + mx + nx 2 )dτ and substituting A = p(1 + m + n) + q + 1 into (3.12), we obtain (still denote τ by t )

where m, p, q, n, a satisfy (3.13). Obviously system (3.14) has the same topological structure as that of system (3.12), since we consider system (3.12) in R + 2 = {(x, y) : x ≥ 0, y ≥ 0}, and 1 + mx + nx 2 > 0 holds for all x ≥ 0. In the following, we study the Hopf bifurcation around Ẽ 2 (1, 1) in system (3.14), which corresponds to the Hopf bifurcation around E 2 (x 2 , y 2 ) in system (2.3).

When conditions in (3.13) hold, system (3.14) has an equilibrium at Ẽ 2 (1, 1). Moreover, (I) when a < a * , Ẽ 2 (1, 1) is an unstable hyperbolic node or focus; (II) when a > a * , Ẽ 2 (1, 1) is a locally asymptotically stable hyperbolic node or focus; (III) when a = a * , Ẽ 2 (1, 1) is a fine focus or center.

Proof. The Jacobian matrix of system (3.14) at Ẽ 2 (1, 1) is

.

Then the determinant of J (Ẽ 2 (1, 1)) is det(J (Ẽ 2 )) = a(m + n + 1)(np + q + 1 − p), and the trace of J (Ẽ 2 (1, 1) ) is

By conditions in (3.13), we can see that det(J (Ẽ 2 )) > 0 and tr(J (Ẽ 2 )) = 0 (> 0 or < 0) if a = a * (a < a * or a > a * ), so the conclusions follow. 2

Next we continue to consider the case (III) of Theorem 3.1 and study Hopf bifurcation around Ẽ 2 (1, 1) in system (3.14). We have the following necessary conditions for Hopf bifurcation:

Firstly we can check the transversality condition d da (trJ (Ẽ 2 )) | a=a * = −(m + n + 1) < 0.

We investigate the nondegenerate condition and stability of the bifurcating periodic orbit from the positive equilibrium Ẽ 2 (1, 1) of system (3.14) by calculating the first Lyapunov coefficient. Let X = x − 1, Y = y − 1, and a = a * . Then system (3.14) can be written as (still denote X, Y, τ by x, y, t , respectively)

where b = p −1 −np. Let ω = b q −b 2 , and make a transformation of x = X, y = 1 q (bX −ωY ) and dt = 1 ω dτ , then system (3.16) becomes (we still denote X, Y by x, y, respectively) 

Using the formula in [33] and calculating first Lyapunov coefficient with the aid of MATLAB, we obtain

where c 0 , c 1 , c 2 , d 0 , d 1 and d 2 are given in (3.9). By conditions in (3.15), the sign of σ 1 is the same as that of

Now we investigate the sign of σ 11 (i.e., σ 1 ) and will see that there exist some values of parameters such that σ 11 = 0 (i.e., σ 1 = 0) when conditions (3.15) are satisfied. We first give the following Lemma. Proof. By conditions (3.15) , it is easy to show that c i > 0 (i = 0, 1, 2) and

Then we have

where m 1 and m 2 are given in (3.9).

(I) Firstly, we consider the discriminant of c 0 + c 1 m + c 2 m 2 = 0, which is given by

where 0 = n 2 > 0, 1 = 8 − 18n + 6n 2 + 4n 3 , 2 = −15 + 60n − 69n 2 + 18n 3 + 6n 4 , 3 = 6 − 48n + 104n 2 − 84n 3 + 18n 4 + 4n 5 , 4 = 1 + 6n − 33n 2 + 52n 3 − 33n 4 + 6n 5 + n 6 .

We consider the third derivative of (p) with respect to p, which is given by (3) (p) = 6 3 + 24 4 p. 3 (1 + n) > 0, we can show that 4 > 0 when 0 < n < 1. Thus

On the other hand, we have p ( 1 1−n ) = 18(n − 1) 2 > 0, and so

Similarly, we have p ( 1 1−n ) = 0 and ( 1 1−n ) = 0, which yield

(II) Secondly, we analyze the sign of m 1 − m 2 , which is same as

By noticing 0 < n < 1 and using a simple analysis, we have

Thus, we can see that m 12 The first derivative of m s1 with respect to p is

Since 0 < n < 1 and p > 1 1−n , one can prove that m s1 (p) < 0. Moreover, we have m s1 ( 1 1−n ) = 0. Thus we have m s1 (p) < 0, i.e., m 1 > − c 1 2c 2 , when conditions (3.15) are satisfied. The sign of − c 1 2c 2 − m 2 is the same as that of

Since 0 < n < 1, by a direct computation, we can show that

It then follows that m s2 (p) < 0, i.e., m 2 > − c 1 2c 2 , when conditions (3.15) are satisfied. 

By a direct computation with conditions in (3.15), we have

where

We can see that c 1 * is monotonic decreasing for 0 < n < 1. Moreover, since c 1 * (0) = 7 and c 1 * (1) = 0, we have c 1 * > 0 (i.e., c m 1 (p * ) > 0) when 0 < n < 1, and thus c(m 1 )(p) > 0, i.e., c(m 1 ) > 0, when p > p * . When 1 1−n < p < p * , we have m 1 < m 2 by the step (II), then m > max{m 1 , m 2 } = m 2 and c(m) > c(m 2 ) by the step (III). Moreover, we can show that c(m 2 ) > 0. In fact, the sign of c(m 2 ) is same as

By direct computation, we have c m2 ( 1 1−n ) = 

We then show that d 2 1 − 4d 0 d 2 = n 2 + 2(n 3 + 3n 2 − 6n + 2)p + n(n 3 + 6n 2 − 15n + 8)p 2 > 0.

In fact, let h(p) = n 2 + 2(n 3 + 3n 2 − 6n + 2)p + n(n 3 + 6n 2 − 15n + 8)p 2 .

By the conditions in (3.15), we obtain

and

Then when 1 3 < n < 1, we have (I) If m ≤ m 3 or m ≥ m 4 (i.e., σ 11 > 0 or σ 1 > 0), then system (3.14) exhibits subcritical Hopf bifurcation and an unstable limit cycle appears around Ẽ 2 (1, 1). (II) If m 3 < m < m 4 and (II.1) q > q 2 (i.e., σ 11 < 0 or σ 1 < 0), then system (3.14) exhibits supercritical Hopf bifurcation and a stable limit cycle appears around Ẽ 2 (1, 1); (II.2) q < q 2 (i.e., σ 11 > 0 or σ 1 > 0), then system (3.14) exhibits subcritical Hopf bifurcation and an unstable limit cycle appears around Ẽ 2 (1, 1); (II.3) q = q 2 (i.e., σ 11 = 0 or σ 1 = 0), then system (3.14) may exhibits degenerate Hopf bifurcation and multiple limit cycles may appears around Ẽ 2 (1, 1).

Here, m 3 , m 4 and q 2 are given in (3.9).

One stable (or unstable) limit cycle arising from supercritical (or subcritical) Hopf bifurcation around the equilibrium Ẽ 2 (1, 1) of system (3.14) is given in Fig. 3.2. In Fig. 3.2(a) , we fix n = 0.1, p = 1.2, m = −0.6 and q = 2, and get a = 0.16 from tr(J (Ẽ 2 )) = 0, and further obtain σ 1 = −3.24949. Finally, we perturb a such that a decreases to 0.12, then Ẽ 2 (1, 1) becomes an unstable hyperbolic focus, leading to a stable limit cycle to appear around Ẽ 2 (1, 1). In Fig. 3.2(b) , we fix n = 0.1, p = 1.2, m = −0.6 and q = 0.6, and get a = 0.16 from tr(J (Ẽ 2 )) = 0, and further get σ 1 = 2.62886. Finally, we perturb a such that a increases to 0.18, and so Ẽ 2 (1, 1) becomes a stable hyperbolic focus, yielding an unstable limit cycle to appear around Ẽ 2 (1, 1). Fig. 3.2(a) , the model has three equilibria, a disease-free equilibrium which is a stable hyperbolic node, two endemic equilibria (one is a hyperbolic saddle and the other is an unstable focus) and a stable limit cycle. If the initial population lies on the right side of the two stable manifolds of the saddle, then the disease will tend to a periodic coexistent oscillation. If the initial population lies on the left side of the two stable manifolds of the saddle, then the disease will die out. Fig. 3.2(b) , the model has three equilibria, a disease-free equilibrium which is a stable hyperbolic node, two endemic equilibria (one is a hyperbolic saddle and the other is a stable focus) and an unstable limit cycle. If the initial population lies on the limit cycle, then the disease will persist in the form of a periodic coexistent oscillation. If the initial population lies inside the limit cycle, then the disease will tend to a stable steady state. If the initial population lies outside the limit cycle, the disease will die out for almost all positive initial populations.

From (II.3) of Theorem 3.5, we know that system (3.14) may exhibit degenerate Hopf bifurcation (i.e., Hopf bifurcation of codimension two) when q = q 2 , m 3 < m < m 4 and conditions (3.15) are satisfied. Using the formal series method in [33] and MATLAB software, we obtain the second Lyapunov coefficient

, (3.20) where

Next we investigate the sign of σ 2 and first give the following lemma. 

By a direct calculation, we have h(m 1 )(p) = 2n + 4(n 2 − 1)p + 3(n − 1) 2 (n + 1)p 2 .

It follows that h(m 1 )(p) > 0 when p > p * . Further, by a direct calculation, we have

so h(m 1 )(p) > 0 (i.e., h(m 1 ) > 0) when p > p * . When 1 1−n < p < p * , we have m 1 < m 2 by the step (II) in the proof of Lemma 3.3. Then m > max{m 1 , m 2 } = m 2 and h(m) > h(m 2 ). Moreover, we can show that h(m 2 ) > 0. In fact

By a direct calculation, we have

It follows that h(m 2 )( 1 1−n ) > 0 and h(m 2 )(p * ) > 0 for 0 < n < 1. Thus, h(m 2 ) > 0 when 1 1−n < p < p * .

Summarizing the above results, we have shown that h(m) > 0 (i.e., 2 +4np +m(np +p +1) > 0) when conditions (3.15) are satisfied. 2

When m > −2 √ n, 0 < n < 1, p > 1 1−n , m 3 < m < m 4 , and by Lemma 3.6, we know that the sign of σ 2 in (3.20) is the same as that of e 0 + e 1 m + e 2 m 2 + e 3 m 3 σ 21 .

(3.21)

Based on extensive numerical calculations, we conjecture that Ẽ 2 (1, 1) is an unstable multiple focus with multiplicity two, i.e., σ 21 > 0 (or σ 2 > 0) when q = q 2 , m 3 < m < m 4 and conditions in (3.15) are satisfied. However, it is very difficult to pursue a pure algebraic proof. Next we combine geometric (graphic) and algebraic methods to show that Ẽ 2 (1, 1) is actually an unstable multiple focus with multiplicity at most two. Actually, we have the following result.

Theorem 3.7. The equilibrium Ẽ 2 (1, 1) in system (3.14) (i.e., E 2 (x 2 , y 2 ) in system (2.3)) is an unstable multiple focus with multiplicity at most two, system (2.3) can exhibit degenerate Hopf bifurcation of codimension two around E 2 (x 2 , y 2 ), and the outer bifurcating limit cycles is unstable.

Proof. The Jacobian of the linearized system of (3.14) shows that Hopf bifurcation occurs at the critical point:

Setting a = a H yields the eigenvalues at the critical point as λ ± = ±ω c i, where

In order to have ω c > 0, it needs 0 < p(1 − n) − 1 < q, (3.22) which in turn yields 1 + m + n > 0 from a H > 0. It is obvious that we need the following condition:

To simplify the analysis, we introduce two new parameters:

Thus, the conditions 0 < p(1 − n) − 1 < q are equivalent to P > 0, Q > 0. Then, a H = P 1+m+n > 0 and ω c is reduced to ω c = √ P Q. Note that the condition q < p a = p a H becomes

Therefore, the conditions in (3.15 ) are equivalent to the following conditions:

Now, introducing the transformation,

and the time rescaling τ 1 = ω c τ into (3.14), we obtain

Next, applying the Maple program [31] we obtain the following focus values:

where F 3 is a lengthy polynomial, which is omitted here for brevity. In the following, we prove that when v 1 = 0, v 2 > 0 for the parameters satisfying the conditions in (3.25) . Since v 1 and v 2 have the same sign of F 1 and F 2 , respectively, we consider the polynomials F 1 and F 2 .

First, suppose that the two polynomial equations F 1 = F 2 = 0 have solutions for P and Q and we study if the solutions satisfy the conditions in (3.25) . Note that F 1 is linear in Q, we solve Q from F 1 = 0 to obtain Q = (1 + m + n)(1 + P ) P (mn + m + 4n) + 2(1 + m + n) (n − 1) P (mn 2 + m 2 + 3mn + 6n 2 − 2n) + m 2 + 4mn + 3n 2 + 2n − 1 , (3.27) and then substitute it into F 2 to get F 2 = 6(n − 1)(1 + m + n) 4 (1 + P ) 3 P P (mn + m + 4n)

Note that the solutions of F 2 = 0 come from the quadratic polynomial equation A 0 P 2 − A 1 P + A 2 = 0. Solving this polynomial equation yields two solutions:

and then we have

To have real solutions for P ± and Q ± , it requires > 0. Summarizing the above results shows that if the two sets of solutions satisfy v 1 = v 2 = 0, they must satisfy the conditions given in (3.25) as well as > 0. Define

It is clear that the conditions given in (3.25) must be in ± . We will show that the obtained solutions P ± and Q ± do not belong to ± , then P ± and Q ± do not satisfy conditions in (3.25) , thus there are no feasible parameter values for v 1 = v 2 = 0, and the best possibility is v 1 = 0 but v 2 = 0, which implies that the maximal number of limit cycles bifurcating from the Hopf critical point is two. To prove the result by combining geometric (graphic) and algebraic methods, first we consider the set − . Define the functions:

then plot the curves = 0 and C k = 0, k = 1, 2, 3 on the m-n plane, and identify the regions on the boundary regime: 0 < n < 1, 1 + m + n > 0, satisfying > 0 and C k > 0, k = 1, 2, 3. If these regions do not overlap, then no solutions exist; if these regions do overlap, we need to further check the last condition in (3.25) , Q < Q max .

For the set − , we show the curves = 0 and C k = 0, k = 1, 2, 3, in Fig. 3 .3(a) with the colors in red, blue, brown and green, respectively. Fig. 3.3(b) is a zoomed area near the m-axis. The region for < 0 is inside the red loop, and > 0 is outside the loop. Thus, we only need to consider the region outside the red loop. We use + or − to indicate C k > 0 or C k < 0. Fig. 3.3(c) is a zoomed area near the corner of (m, n) = (−2, 1), showing that all three curves are tangent to the line 1 + m + n = 0 at the corner. This can be rigorously proved using an algebraic argument. Figs. 3.3(a), (b) and (c) clearly show that there are no regions with > 0 and all three functions C k > 0, k = 1, 2, 3. Thus, the set − does not have feasible solutions satisfying v 1 = v 2 = 0.

Next, consider the set + . We have a similar result, as shown in Fig. 3 .4. Again, this set also does not have feasible solutions satisfying v 1 = v 2 = 0.

To this end, we have shown that under the restrictions on the parameters, as given in (3.25), it is not possible to have solutions such that v 1 = v 2 = 0 and to have three limit cycles bifurcating from the equilibrium (x 0 , y 0 ) = (1, 1) due to Hopf bifurcation. Therefore, the next best choice is to have v 1 = 0 and v 2 = 0, and P can then be treated as a free parameter. The above proof implies that v 2 must be kept the same sign for the feasible parameter values. Otherwise, if it can change sign, then it must have solutions for v 2 = 0, which contradicts the above result. To determine whether the sign of v 2 is positive or negative, we only need to choose a solution such that v 1 = 0 and then calculate v 2 . For example, taking P = n = m = 1 10 , we obtain

giving v 1 = 0 and > 0, implying that v 2 > 0 for any parameter values which satisfy the conditions given in (3.25) . Thus, the outer bifurcating limit cycle is always unstable and the inner one is stable, both enclose an unstable focus (1, 1) . The proof is complete. 2 Remark 3.3. If m = 0 (i.e., β = 0), then σ 21 = e 0 > 0 (i.e., σ 2 > 0), thus Ẽ 2 (1, 1) (i.e., E 2 (x 2 , y 2 )) is an unstable multiple focus with multiplicity two when q = q 2 , m 3 < m < m 4 and conditions in (3.15) are satisfied. Moreover, when β = 0 (i.e., m = 0), our system (2.2) is reduced to system (1.2) in [25] , where the authors only discussed the Hopf bifurcation of codimension one, thus our results about Hopf bifurcation of codimension two can be seen as a complement of [25] by considering a more general system. On the other hand, when β = 0 (i.e., m = 0), our system (2.2) also becomes system (1.3) in [27] , where the authors showed the existence of degenerate Hopf bifurcation by complex algebraic calculations, our Theorems 3.5 and 3.7 give more explicit conditions and our proofs are easier to follow.

Next we give some numerical simulations in Fig. 3 .5 to show the existence of two limit cycles based on Theorems 3.5 and 3.7. Firstly, we fix n = 1/3, p = 2 and m = −1, then get q = 23 15 from σ 1 = 0, and get a = 1 from tr(J (Ẽ 2 )) = 0, finally get σ 2 = 2.65719, i.e., Ẽ 2 (1, 1) is an unstable multiple focus with multiplicity two for those fixed parameters. Next we first perturb q such that q increases to 23 15 + 0.1, then Ẽ 2 (1, 1) becomes a stable multiple focus with multiplicity one, an unstable limit cycle occurs around Ẽ 2 (1, 1) which is the outer limit cycle in Fig. 3 .5. Secondly, we perturb a such that a decreases to 1 − 0.003, then Ẽ 2 (1, 1) becomes an unstable hyperbolic focus, another stable limit cycle occurs around Ẽ 2 (1, 1), which is the inner limit cycle in Fig. 3 .5. Fig. 3 .5, we can see that the existence of multiple periodic coexistent oscillations and coexistent steady states when the psychological effect α is smaller than the critical value α 0 but the infection rate k is greater than the first critical value k 0 , i.e., n < n * and A > mp. The disease will die out for almost all initial populations outside the outer unstable periodic orbit, and will tend to periodic outbreaks for almost all initial populations on or inside the outer unstable periodic orbit, and will persist in the form of positive steady states when the initial population lies on the positive equilibria or their stable manifolds.

In section 1 we mentioned that nonmonotone and saturated incidence rates may apply to some infectious diseases such as influenza. In this section, we use model (2.1) to simulate the reported human influenza data from China. From the Chinese Center for Disease Control and Prevention (CCDC) [8] , we obtain the annual data on human influenza cases from 2004 to 2017 which are shown in Table 1 . Most parameter values can be obtained from the literature or by estimation. We estimate k, α, β, δ and R(0) by using the unconstrained optimization functions f minsearch, a part of the optimization toolbox in MATLAB, to fit I (t i ) through discretizing the ordinary differential system (2.1) as follows

kS(t i )I (t i ) 2 1 + βI (t i ) + αI (t i ) 2 − (d + μ)I (t i ) t + I (t i ).

The fitting is to minimize the objective function Table 1 The data on human influenza cases in China from 2004-2017 ( [8] ). 

where I (t i ) is numerical solution of model (2.1) and Î (t i ) is the reported human influenza data from China.

The parameter values are listed in Table 2 . From the data in 2004, we know that I (0) = 49496; we make the data fitting to obtain that R(0) = 13, then S(0) = 1.29983 × 10 9 ; and assume that μ = 365 5 from Casagrandi et al. [7] , data fitting gives δ = 0.33042, k = 2.8668 × 10 −8 , α = 4.9077 × 10 −6 and β = 78.12.

Based on the parameter values given in Table 2 , we use model (2.1) to simulate the data from 2004 to 2017 and predict the trend of human influenza infections in China. Fig. 4 .1 represents the simulation of our model with reasonable parameter values which provides a good match to the data on infected human influenza cases in China from 2004 to 2017. In Fig. 4 .2 our model predicts that the number of human influenza infection will increase in the next a couple of years.

In most classic SIRS epidemic models, the incidence rate is either monotonic or nonmonotonic, and the infection force g(I ) usually tends to zero as I → ∞ when the incidence rate is non-monotonic, which may be unreasonable for some specific infectious diseases. In this paper, we considered the SIRS epidemic model (1.1) with a generalized nonmonotone and saturated incidence function g(I ) = kI 2 1+βI +αI 2 , which exhibits three different properties: (i) when β ≥ 0, g(I ) is monotonic, which always increases and then tends to a saturated level k α as I intends to infinite. (ii) When −2 √ α < β < 0, g(I ) is nonmonotonic, which increases when I is small and decreases when I is large, and finally tends to a saturated level k α as I tends to infinite (see Fig. 1.2) , which can be used to describe the phenomenon that the infection force first increases to a maximum when a new infectious disease emerges, and then decreases due to psychological effect and eventually tends to a saturation level due to crowding effect. (iii) When β > −2 √ α, g(I ) has both properties of (i) and (ii), i.e., g(I ) is first nonmonotonic and then monotonic when β increases from negative to positive. In this paper, we performed qualitative and bifurcation analysis, and found that the model (1.1) with general incidence rate (1.8) has complicated dynamical behaviors and bifurcation phenomena. It has been shown that there exists a weak focus of multiplicity at most two and a cusp of codimension at most two for various parameter values, and the model undergoes saddle-node bifurcation, Bogdanov-Takens bifurcation of codimension two, Hopf bifurcation and degenerate Hopf bifurcation of codimension two as the parameter values vary. By considering a more general model, our results show that Bogdanov-Takens bifurcation of codimension higher than two cannot occur, which clarifies and corrects the corresponding results by Tang et al. [27] . Moreover, our results about degenerate Hopf bifurcation (i.e., Hopf bifurcation of codimension two) can also be seen as a complement of the results in Ruan and Wang [25] . We also used model (2.1) to simulate the reported human influenza data from 2004 to 2017 reported by the Chinese Center for Disease Control and Prevention.

In recent years, people realize the importance of the inhibition effect from the behavioral change of susceptible individuals when their number increases or from the crowding effect of the infective individuals, and the psychological effect or inhibition effect can cause people to adopt some aggressive measures and policies, such as border screening, mask wearing, quarantine, isolation and so on, which have been proved to be very effective in reducing the contact rate or the infection force. However, it is difficult to measure the psychological effect. So can we give a quantitative description about the psychological effect with respect to variation of the parameters involved in the infectious diseases? In this paper, we presented a concrete critical value to measure the psychological effect. More specifically, we have shown that there exists a critical value α = α 0 for the psychological or inhibitory effect and two critical values k = k 0 , k 1 (k 0 < k 1 ) for the infection rate such that: (i) when α > α 0 , or α ≤ α 0 and k ≤ k 0 , the disease will die out for all positive initial populations; (ii) when α = α 0 and k 0 < k ≤ k 1 , the disease will die out for almost all positive initial populations; (iii) when α = α 0 and k > k 1 , the disease will persist in the form of a positive coexistent steady state for some positive initial populations; and (iv) when α < α 0 and k > k 0 , the disease will persist in the form of multiple positive periodic coexistent oscillations and coexistent steady states for some positive initial populations.

Periodicity in an epidemic model with a generalized nonlinear incidence

Bifurcation analysis of an SIRS epidemic model with generalized incidence

Bifurcations of a limit cycle for a family of vector fields on the plane

Versal deformations of a singular point on the plane in the case of zero eigen-values

Serio, I modelli matematici nella indagine epidemiologica. I Applicazione all'epidemia di colera verificatasi in Bari nel 1973

A generalization of the Kermack-McKendrick deterministic epidemic model

The SIRC model and influenza A

Chinese Center for Disease Control and Prevention, Reports on Infectious Diseases

Homoclinic orbits in a disease transmission model with nonlinear incidence and nonconstant population

A disease transmission in a nonconstant population

Global dynamics for a predator-prey system of Leslie type with generalized Holling type III functional response

The mathematics of infectious disease

Some epidemiological models with nonlinear incidence

Bifurcations of an SIRS epidemic model with nonlinear incidence rate

Bifurcation analysis in a predator-prey model with constant-yield predator harvesting

A contribution to the mathematical theory of epidemics

Bifurcation of an SIS model with nonlinear contact rate

Global dynamics of avian influenza epidemic models with psychological effects

Dynamical behavior of epidemiological models with nonlinear incidence rates

Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models

Multiparametric bifurcations for a model in epidemiology

Bifurcations of an epidemic model with non-linear incidence and infectiondependent removal rate

Differential Equations and Dynamical Systems

Dynamical behavior of an epidemical model with a nonlinear incidence rate

Forced oscillations and bifurcation, in: Applications of Global Analysis I

Coexistence of limit cycles and homoclinic loops in a SIRS model with a nonlinear incidence rate

Epidemic models with nonlinear infection forces

Global analysis of an epidemic model with a nonlinear incidence rate

Qualitative analysis of an epidemic model

Computation of normal forms via a perturbation technique

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

Qualitative Theory of Differential Equations