BIFURCATION ANALYSIS OF A DISCRETE SIS MODEL WITH BILINEAR INCIDENCE DEPENDING ON NEW INFECTION

A discrete SIS epidemic model with the bilinear incidence depending on the new infection is formulated and studied. The condition for the global stability of the disease free equilibrium is obtained. The existence of the endemic equilibrium and its stability are investigated. More attention is paid to the existence of the saddle-node bifurcation, the flip bifurcation, and the Hopf bifurcation. Sufficient conditions for those bifurcations have been obtained. Numerical simulations are conducted to demonstrate our theoretical results and the complexity of the model.

1. Introduction.Differential equations and difference equations are widely applied in epidemiological modeling.They are two typical mathematical approaches for modeling infectious diseases.Since the theory and method for dynamical studies of differential equations have developed much more completely than those for difference equations, there are relatively less difference equations in epidemiological modeling compared with differential equations.In recent years, there were increasing interest and research results on discrete epidemic models [1, 2, 10-12, 14, 30, 32].The fact that the epidemiological data are usually collected in discrete time units, such as daily, weekly or monthly, makes the discrete model a natural choice to describe a disease transmission.The straightforward recurrence relationship of the difference equation models is easier to be understood, which is also a prominent advantage over the differential equation models.The direct comparison of the model results with the actual data provides us a fast and simple way to validate the model structure and parameter estimation.The fact that the discrete models exhibit richer dynamical behavior than the continuous models brings more challengeable problems for researches, and more interesting results can be obtained.For example, the simple logistic model x n+1 = rx n (1 − x n /K), Ricker model x n+1 = x n e r−xn/K [23][24][25], and Hassell model x n+1 = λx n (1 + ax n ) −b [18] exhibit rich dynamical behavior.

HUI CAO, YICANG ZHOU AND ZHIEN MA
There are increasing interest and more studies on the discrete time epidemic models recently.Various discrete epidemic models have been successfully applied to describe the infectious disease transmission, such as SARS, tuberculosis, HIV/AIDS [6-8, 28, 29, 31].The theoretical study of discrete epidemic models focus on the computation of the basic reproduction number [3,9,22], the existence and the global stability of the disease free equilibrium [7,8,11,12,28], the existence and local stability of the endemic equilibrium [4,15], and the persistence of the disease [7,8].The attention has also be paid to various bifurcations of the discrete epidemic models, the equilibrium bifurcation [10,[19][20][21][22], the transcritical bifurcation, flip bifurcation, saddle-node bifurcation, Hopf bifurcation, and the bifurcation to chaos [10][11][12].
Similar to the continuous epidemic models, most of the discrete epidemic models have the compartment structures.The total population is grouped into different epidemiological compartments, such as the susceptibles, the infectives, and the removed.The key step in formulating discrete epidemic models is to find the recurrent relationships among those epidemiological compartments between time t and time t + 1.Two main approaches are frequently used to describe the population movement among those compartments.The first one is analogous to the discretization of the corresponding continuous models [1,2], and the second one is to seek the probability from one compartment to another [10][11][12].The transmission rate, the death rate, and recovery rate used in the first approach make those kinds of models can be directly applied to a specific infectious disease.The probabilities used in the second approach make the solutions of those models often positive.
Following the modeling approach by Allen and Burgin [1,2], we propose a discrete SIS model and study its dynamical behavior.In Section 2, we present the SIS model with the bilinear incidence depending on the new infection and simplify it.We define the basic reproduction number and the invariant domain for the model.In Section 3, we discuss the existence and the number of the endemic equilibria.In Section 4, we study the global stability of disease free equilibrium, and the local stability of the endemic equilibria.In Section 5, we investigate the flip bifurcation, saddlenode bifurcation and Hopf bifurcation of the model.Numerical simulations are conducted to demonstrate our theoretical results and show the complexity of the model dynamics.Concluding remarks and discussions are given in the last section.
2. The SIS model and its basic reproductive number.The incidence plays a crucial role in dynamics of epidemic models.The bilinear incidence and the standard incidences are two simple and frequently used ones in continuous and discrete models.We consider a discrete epidemic model with SIS structure and use the bilinear incidence depending on the new infection.The bilinear incidence depending on the new infection is a better choice to describe the influence of control measures.Let N (t) represent the number of the total population in a community at time t.We divide the total population into susceptible and infectious compartments, where S(t) and I(t) denote the numbers of individuals in the susceptible and infectious compartments at time t, respectively.The discrete SIS model with bilinear incidence depending on the new infection is S(t + 1) = S(t) + Λ − βS(t)I(t)e −kI(t−1)S(t−1) − µS(t) + γI(t), I(t + 1) = I(t) + βS(t)I(t)e −kI(t−1)S(t−1) − (µ + γ)I(t), where Λ is the constant recruitment of the population, µ is the per capita natural death rate, and γ is the per capita recovery rate.All parameters are positive, and 0 < β < 1 and 0 < µ + γ ≤ 1.The number of new infection between time interval [t, t + 1] is βS(t)I(t)e −kI(t−1)S(t−1) .The factor e −kI(t−1)S(t−1) describes the fact that effective control measures will be implemented to reduce the transmission of the disease.The impact of the control measures depends on the new infection: more infection will lead to more stringent control, which will reduce the infection.Following the next generation matrix approach [3], we know that the basic reproductive number, R 0 , of model ( 1) is R 0 is the average number of secondary cases produced by a typical single infected case in a completely susceptible population.R 0 is a useful quantity to predict whether or not the disease can spread: if R 0 < 1, the infection will ultimately die out, while if R 0 > 1, the infection will keep persistence in the population.If R 0 is large, it will be very harder to control the disease spread and a major epidemic may outbreak.
Let N (t) = S(t) + I(t), adding all equations in (1), we can obtain the equation for the total population N (t) It is clear that equation ( 2) has a unique equilibrium N * = Λ µ , which is globally asymptotically stable.That is, lim the second equation of model (1) leads to the following limiting system I(t + 1) = I(t) + β N * − I(t) I(t)e −kI(t−1) N * −I(t−1) − (µ + γ)I(t). (3) The limiting system of model ( 1) is one dimensional difference equation with delay.
After defining x 1 (t) = I(t), y 1 (t) = I(t − 1), and N = N * , model (3) can be rewrote as The qualitative equivalence of a dynamical system and its limiting system has been established by Thieme [26,27] for continuous systems.We study the dynamical behavior of model (4) since model (4) exhibits the same qualitative dynamics as those of system (1).From the expression of model (2) and model (3), it is not difficult to obtain the dynamical equivalence of those models.After introducing new variables x(t) = x1(t) N , y(t) = y1(t) N , and parameters a = βN , b = kN 2 , c = µ + γ, model (4) can be simplified to The simplified model (5) has less parameters and is easier to analyse.The epidemiological interpretation requires that 0 < c < 1 and 0 ≤ x(t) ≤ 1.A sufficient condition for 0 ≤ x(t) ≤ 1 is that the parameters c and a locate in the domain The condition (c, a) ∈ P + is very strong to ensure that the solution of (5) satisfies 0 ≤ x(t) ≤ 1 and 0 ≤ y(t) ≤ 1 if 0 ≤ x(0) ≤ 1 and 0 ≤ y(0) ≤ 1.We assume that (c, a) ∈ P + when we discuss the stability of equilibrium of model (5).The requirement (c, a) ∈ P + may be released when we investigate the bifurcation of model (5).In the remaining part of this paper, we will focus on the study of model (5).The stability and bifurcation results of model (5) will equivalently hold true for the original model.

3.
The existence of equilibria.This section deals with the existence of the equilibria of model (5).It is clear that E 0 (0, 0) is an equilibrium point of model (5).There also exists the positive equilibrium points E(x, y) when c < a, which implies that R 0 > 1. Namely, when R 0 > 1, model (5) has the positive equilibrium points E(x, y), and the positive equilibrium points E(x, y) satisfies the following equations: We denote Obviously, f (x) ≥ 0 for x ∈ [0, 1], f (0) = a, and f (1) = 0. Furthermore, we have .
From the fact that we obtain the following results: 2 ) .In the case where R 0 > 1 (which is equivalent to c < a), the existence and the number of the positive equilibria of model ( 5) are given in the following theorem.Theorem 3.1.When R 0 > 1 and b − 8 ≤ 0, model (5) has only one positive equilibrium point E 1 (x 1 , x 1 ).When R 0 > 1 and b − 8 > 0, there are these cases: 3 ), and Theorem 3.1 can be proved directly by using the intermediate value theorem for continuous function.We omit the detailed process and demonstrate the idea by  5) for different parameters 4. The stability of equilibria.Let E(x, x) be a positive equilibrium point of model (5).The linearization matrix of (5) at the equilibrium point E(x, x) is The characteristic equation of matrix A is Let λ 1 and λ 2 be the two roots of Eq. ( 6), we use the following definitions and conclusions [19,20]. ( In the following, we will discuss the stability of the disease free equilibrium and endemic equilibria of model ( 5).
Proof.The condition R 0 < 1 implies that 0 < a < c < 1.When E(x, x) = E 0 (0, 0), we have In this case, the eigenvalues of matrix A are λ 1 = 1−c+a and λ 2 = 0, respectively.It is obvious to see that Therefore, E 0 is locally asymptotically stable when R 0 < 1.
In order to prove the global stability of the disease free equilibrium we rewrite model (5) to be From Eq. ( 7) and the fact 0 ≤ x(t) ≤ 1, we have The variable substitution z(t) = ax(t) 1 − c + a leads to the comparison equation The famous result on the dynamics of the discrete population model implies that ) has a zero equilibrium point, which is globally asymptotically stable if 0 ≤ r < 1.Therefore, the comparison theorem implies that the zero solution of Eq. ( 7) is globally asymptotically stable if 0 < 1 − c + a < 1.From the equivalence of 0 < 1 − c + a < 1 and R 0 < 1, we know that the disease free equilibrium E 0 of model ( 5) is globally asymptotically stable when R 0 < 1.

4.2.
The stability of the endemic equilibria.Next, we discuss the stability of the endemic equilibria of model (5).The equation a(1 − x)e −bx(1−x) = c of the positive equilibrium can simplify the linearized matrix A of model ( 5) at the equilibrium E(x, x) The characteristic equation of matrix A is Following results give the relationship between roots and coefficients of a quadratic equation.(5) λ 1 and λ 2 are the conjugate complex roots with The function g 2 (x) is a cubic polynomial and satisfies is the unique positive root of g 2 (x) = 0 in the interval (0, 1).From the continuity of g 2 (x) we know that g 2 (x) > 0 if x ∈ (0, x 2 g ), g 2 (x) = 0 if x = x 2 g , and g 2 (x) < 0 if x ∈ (x 2 g , 1).Consequently, h(−1) > 0 if x ∈ (0, x 2 g ), h(−1) = 0 if x = x 2 g , and h(−1) < 0 if x ∈ (x 2 g , 1).Based on the analysis of h(1) and h(−1), we can prove following stability results.Theorem 4.3.If R 0 > 1 and b < 8, then the stability results of the unique endemic equilibrium E 1 (x 1 , x 1 ) of model (5) are as follows: (i) E 1 is a sink and is locally asymptotical stable if The function −bc(2x 1 − 1)x 1 has its maximum bc 8 at x 1 = 1 4 .Consequently, −bc(2x 1 − 1)x 1 < 1 for x 1 ∈ (0, 1).From the aforementioned results on h(−1) and Lemma 4.2, we know that the conclusions given in Theorem 4.3 hold true.
When b = 8, the similar stability result of E 1 (x 1 , x 1 ) can be obtained if More careful analysis should be done for the case x 1 = 3 4 .
When b > 8 more endemic equilibria may appear and the stability analysis will be more complicated.We use following special case to show the number and stability of endemic equilibria of model (5).When we take b = 9 and a = 1, then x * 1 = 0.6667, x * 2 = 0.8333, c 1 = 0.04511, and c 2 = 0.04775.The existence and stability of endemic equilibria are given in the following statements.
Assume that b = 9 and a = 1, then (i) Model ( 5) has one endemic equilibrium is a sink and it is locally stable.(iii) Model ( 5) has three endemic equilibria E 1  3 , E 2 3 , and 3 is a saddle, and E 3 3 is a sink and it is locally stable.
These statements are obtained numerically.In fact, we have tried different parameter values of a and b to investigate the number and stability of the endemic equilibria of model (5).We have got similar results.
5. The bifurcation.Bifurcation may lead to different dynamical behaviors of a model when parameters pass through a critical value.Bifurcation usually occurs when the stability of an equilibrium changes.In this section, we will discuss the flip bifurcation, the saddle-node bifurcation, and the Hopf bifurcation of model (5).
5.1.The flip bifurcation.We investigate the flip bifurcation under the condition b < 8 and R 0 > 1.In this case, model (5) has the unique endemic equilibrium E 1 (x 1 , x 1 ).The stability analysis in Section 4 shows that the endemic E 1 has a eigenvalue −1 when x 1 = x 2 g , and E 1 (x 2 g , x 2 g ) is non-hyperbolic.The flip bifurcation can occur in the neighborhood of the endemic equilibrium E 1 (x 2 g , x 2 g ) when parameters pass through a critical point.
The calculation of flip bifurcation is quite complicated and we take a = 3ce 2b 9 to make the process simple.In the case, the endemic equilibrium 3 .The linearization matrix of (5) at the equilibrium point and the characteristic equation of matrix A is )( 13 ỹ(t)+ỹ 2 (t)) − 2 3 c, μ(t + 1) = μ(t), ỹ(t + 1) = x(t).( 9) Taylor expansion of model ( 9) at (x, ỹ, μ) = (0, 0, 0) is where Let Using the transformation for model (10) yields where From the center manifold theory of discrete system we know that there exists a local manifold of model ( 12) [5].The local manifold has the following expansion After substituting the expansion into model ( 12) and using the invariant property of the local manifold, the straightforward and careful calculation gives From the second equation of model (12) we know that µ(t) is always constant.Therefore, the one dimensional model induced by the center manifold is where It is not difficult to verify that G(0, µ)=0, ∂G(0, 0) . Numerical computation shows that ∂ 3 G 2 (0, 0) ∂u 3 < 0 when c ∈ (0.782, 1) (see Fig 2(a)).Therefore, model (5) will undergo a flip bifurcation at E 1 2 3 , 2 3 and the bifurcation solution of period 2 is stable [16].
The period 2 solution from the flip bifurcation of model ( 5) is shown in Figure 2 The saddle-node bifurcation, or the fold bifurcation, describes the phenomenon that a non-hyperbolic equilibrium separates into two equilibria, one is a saddle, and the other is a node.In this subsection, we will prove the existence of the saddle-node bifurcation by a direct calculation of the number of the equilibria and their stability though there exist the standard theory and schedule to determine the saddle-node bifurcation.
From the existence result in Section 3, we know that ) is an positive equilibrium of model ( 5) if b > 8 and c = c 1 .The stability analysis in Section 4 shows that one eigenvalue of the linearization matrix of (5) at the equilibrium ) is 1, which implies that the model ( 5) may undergo the saddle-node bifurcation.Following theorem conforms the saddle-node bifurcation at the equilibrium E * 1 (x * 1 , x * 1 ).Theorem 5.2.If b > 8, then model (5) will undergo a saddle-node bifurcation at when the parameter c increases from the critical value c 1 , where Proof.When b > 8 and = c 1 , we see that ) is an positive equilibrium of model (5).Further calculation shows that the eigenvalues of the linearized matrix of model ( 5) at E * 1 (x * 1 , x * 1 ) are Let u be a small positive number.The straightforward calculation shows that x * 1 − u and x * 1 + u are two roots of the equation f (x) = c, where, From relationship between the equilibrium of model ( 5) and the root of the equation are the positive equilibria of model (5).The equilibria E 1  3 and E 2 3 are very close to E * 1 when u is sufficiently small.Those two equilibria will collide and become E * 1 (x * 1 , x * 1 ) when u tends to zero.
Next, we will determine the stability and the type of those two positive equilibria E 1  3 and E 2 3 .We give the detailed analysis for E 2  3 , and the analysis of E 1 3 is omitted since the process is similar.The characteristic equation of the linearized matrix at Let λ 1 (u) and λ 2 (u) be two roots of equation ϕ(u, λ) = 0. From the expression ϕ(u, λ) we know that the eigenvalue, λ 1 (u) or λ 2 (u), is continuously differentiable with respect to u, satisfying The direct calculation gives that By using the derivative rule for implicit functions we have ∂c ∂u = 0, and The last inequality holds true since , and Therefore, λ 1 (u) > 1 if u > 0 and u is sufficiently small.0 < λ 2 (0) < 1 and dλ du u=0 > 0 say that 0 < λ 2 (u) < 1 for sufficiently small u.The fact that λ 1 (u) > 1 and 0 < λ 2 (u) < 1 implies that the positive equilibrium E 2 3 (x * 1 + u, x * 1 + u) is an unstable saddle.The similar argument and process show that the positive equilibrium is a stable node for sufficiently small u.Combining those results together we know that saddle-node bifurcation will be bifurcated from the equilibrium E * 1 (x * 1 , x * 1 ) when c increases from c 1 .The proof of Theorem 5.2 completes.
Remark 1.If b > 8, the similar arguments can confirm that (5) will undergo the saddle-node bifurcation at E * 2 (x * 2 , x * 2 ) when the parameter c decreases from the critical value c 2 , where The saddle-node bifurcation of model ( 5) is shown in Figure 3. Parameter values a = 0.008 and b = 9 are same for those two subplots, whereas c = c 1 = 0.036089 for subplot (a) and c = 0.036094 for subplot (b), respectively.Figure 3 (a) shows that E * 1 is a non-hyperbolic equilibrium.The solutions initiating from the left part of the neighborhood tend to E * 1 when time goes to infinity, whereas the solutions initiating from the right part of the neighborhood leave E * 1 when time goes to infinity.Figure 3 (b) shows that E * 1 separates into a stable node E 1 3 (0.6611, 0.6611) and and a saddle E 2 3 (0.6722, 0.6722) when c increases a little from the critical value c 1 .Solutions starting in the neighborhood of E 1 3 tend to it when time goes to infinity, whereas solutions starting in the neighborhood of E 2 3 leave the neighborhood when time goes to infinity.for the discrete models is similar to that of continuous models.In this subsection we discuss the Hopf bifurcation of model ( 5) when b > 8 and a = 4  3 ce 3b 16 .In this case E 1 3 1 4 , 1  4 is a positive equilibrium of ( 5), and the linearization matrix of (5) at the equilibrium point E 1  3 is The characteristic equation of matrix A is The Taylor expansion of model ( 16) at (x, ỹ) = (0, 0) is  4 ). ( The eigenvalues of the linearized matrix are Let b be the bifurcation parameter and b 0 = 8 c .The direct calculation yields that where From Theorem 3.5.3 of [17] we know that the existence of Hopf bifurcation can be determined by the quantity θ, where, By a straightforward and tedious calculation we obtain that Using the Hopf bifurcation theorem in [17] we obtain that there exists a Hopf bifurcation when b passes through b 0 .The proof completes. The solution curve and the invariant curve of Hopf bifurcation are shown in  5) at ( 14 , 1 4 ) is the positive equilibrium point of model ( 5).The characteristic equation of the linearization matrix of (5) at the equilibrium It is easy to have that h(1) = 6.Conclusion and discussion.We have studied a discrete SIS model with the bilinear incidence depending on the new infection.The incidence can make the model more realistic to describe an infectious when the control measures become more stringent if the new infection increases.The global stability of disease free equilibrium, the existence and local stability of endemic equilibrium, the flip bifurcation, the saddle-node bifurcation, and the Hopf bifurcation of the model have been studied.Although the bifurcation results in this paper are given under certain assumptions, numerical simulations shows that those bifurcation hold true for more values of parameters.
The bifurcation results of model (5) indicate that the simple model exhibits very complex dynamical behaviors.Although those bifurcations are investigated separately in the paper, it is very possible that two bifurcations may happen together at one equilibrium, that is, if the eigenvalues of the linearized matrix are 1 and -1, then more complicated bifurcation may appear.Besides those three kinds of bifurcation in our analysis, we can also find the chaotic behavior for certain parameters.The numerical simulation, given in Figure 5, shows that model (5) exhibits chaotic behavior for b < 8. From Figure 5 we see that the endemic equilibrium bifurcates at R 0 = 1.The endemic equilibrium is stable when R 0 is not large enough, then the endemic equilibrium becomes unstable, and a stable periodic solution of period two appears as R 0 increases.The periodic solution loses its stability when R 0 passes certain value and a stable periodic solution of period four bifurcates.The bifurcation process may continue to chaos with the increase of the basic reproductive number R 0 .5) can undergo the forward-backward-forward bifurcation for R 0 > 1. Model (5) has only one endemic equilibrium when 1 < R 0 < R 1 0 , two endemic equilibria when R 0 = R 1 0 , three endemic equilibria when R 1 0 < R 0 < R 2 0 , two endemic equilibria when R 0 = R 2 0 , and one endemic equilibrium when R 0 > R 2 0 (see Figure 6).In the case R 0 = R 1 0 or R 0 = R 2 0 , one of the endemic equilibria is not hyperbolic, and saddle-node bifurcation may occurs when R 0 increases from R 1 0 or decreases from R 2 0 .Bi-stability is found when R 1 0 < R 0 < R 2 0 .Global stability of the endemic equilibrium is a challenging problem in dynamical studies of epidemic model.We have not paid enough attention to the global stability of the endemic equilibrium of model (5).The global stability is only established for

2 − 4 4 .
< 0. By Lemma 4.2, we know that the eigenvalues of the equation h(λ) = 0 are complex, and satisfy that |λ 1 | = |λ 2 | = 1.Therefore, model (5) may undergo the Hopf bifurcation at the endemic equilibria E The possible Hopf bifurcation consists of a surface in parameter space: a =