Bifurcation analysis of a discrete SIR epidemic model with constant recovery

We state and study a discrete SIR epidemic model with bilinear incidence rate and constant recovery. We obtain conditions for the existence of the disease-free equilibrium and endemic equilibria. Theoretic analysis shows that the disease-free equilibrium is globally asymptotically stable when the basic reproduction number is less than unity, and the numerical simulations illustrate that it is asymptotically stable when the number is greater than unity. We also obtain conditions for the stability of the endemic equilibria. More attention is paid to the existence of the saddle-node bifurcation, the flip bifurcation, the $1:1$1:1 resonance, and the Neimark–Sacker bifurcation. We obtain sufficient conditions for those bifurcations. Our numerical simulations demonstrate our theoretical results and the complexity of the model.


Introduction
Mathematical models describing the population dynamics of infectious diseases have been playing an important role in better understanding epidemiological patterns and disease control for a long time [1]. In some cases a simple mathematical model can reveal the nature of the infectious disease transmission, which plays an important role in the control and prevention of the infectious diseases. In general, differential equations and difference equations are two typical mathematical approaches to modeling epidemic dynamical systems.
In recent years, there has been an increasing interest on discrete population dynamical systems because of their rich dynamics behavior and suitability. On the one hand, the discrete models exhibit richer dynamical behavior than the continuous models, which brings more challengeable problems for researches, and more interesting results can be obtained. For example, the simple logistic model x n+1 = rx n (1x n /K), the Ricker model x n+1 = x n e r-x n /K , and the Hassell model x n+1 = λx n (1 + ax n ) -b all exhibit rich dynamical behavior [2][3][4][5]. On the other hand, since epidemiological data are usually collected at discrete time intervals, it is easier to compare the numerical results of discrete epidemic model with real data. In addition, the straightforward recurrence relationship of the dif-ference equation models is easier to be understood, which is also a prominent advantage over the difference equation models.
In theoretical epidemiology the Euler discretization has been extensively utilized to build the discrete epidemic models [1,[6][7][8][9][10][11][12][13][14]. The main idea is discretizing the existing continuous models. However, the disadvantage of this method is that it cannot guarantee the nonnegativity and boundedness of the solutions of the systems in some cases. As a result, some scholars establish discrete epidemic models by using the probability method [15][16][17]. Although the probability method can guarantee the nonnegativity and boundedness of the solutions, the model is more complex, and theoretical analysis is difficult. In this paper, we build a discrete SIR model by applying the Euler discretization. We are lucky that the solution of our model is nonnegative and bounded.
Theoretical results of discrete epidemic models mainly focus on these aspects, including the computation of the basic reproduction number [8,18,19], the comparison between continuous-time epidemic models and discrete-time epidemic models [10,20], the local stability and global stability of the disease-free equilibrium and endemic equilibrium [6,7,9,11,15,16,21], the extinction, persistence, and permanence of a disease [12,[22][23][24][25], periodic systems [8,19,26], bifurcations and chaos phenomena [1,7,9,13,14,17,25,27,28], and so on. In particular, these papers about the bifurcation analysis of discrete epidemic model mainly involve the conditions on the existence of codimension-one bifurcations, such as fold bifurcation, flip bifurcation, Neimark-Sacker bifurcation [1,7,9,25,27,28], and so on, which were derived by using the center manifold theorem and bifurcation theory. Here [7] and [9] both focus on the influence of discretization step size on the dynamic behavior of the system. However, there are few studies on codimension-two bifurcations [13,14], for example, bifurcations associated with 1 : 1 strong resonance. In [14], based on the analysis of various strong resonance, a controller was designed to eliminate the disease. In this paper, we are more interested in finding conditions of the existence of the saddle-node bifurcation, the flip bifurcation, the 1 : 1 resonance, and the Neimark-Sacker bifurcation.
Wang and Ruan [29] have used a continuous-time SIR model tell us that it is unnecessary to take such a large treatment capacity that endemic equilibria disappear to eradicate the disease, and the outcome of disease spread may depend on the position of the initial states for certain range of parameters. We put the continuous-time SIR model in [29] as follows: where S(t), I(t), and R(t) denote the numbers of susceptible, infective, and recovered individuals at time t, respectively, Λ is the recruitment rate, β is the transmission rate, d is the natural death rate, γ is the spontaneous recovery rate of the infective individuals, and h(I) is the removal rate of infective individuals due to the treatment of the form This means that if the number of infective individuals is greater than 0, then it is possible for the infective individuals to be cured, leaving the infective class and entering the recovered class. If the number of infective individuals were 0, then no infective individuals would enter the recovered class.
In this paper, we apply the forward Euler scheme to model (1.1) and obtain the following discrete-time SIR epidemic model: Here all parameters have the same meaning as in system (1.1), and h(I) = m ≥ 0. We assume that all the parameters are positive constants except for m. Based on the epidemiological meaning, we require that 1d > 0, 1dγ > 0, and Λ > m.
Due to the biological interpretation of system (1.3), only nonnegative solutions are meaningful to be considered. The following result reveals that the solutions (S(t), I(t), R(t)) of system (1.3) with non-negative initial value ultimately remain nonnegative and bounded. (1.5) Obviously, systems (1.5) and (1.3) are equivalent. In the following, we show that the solutions (N(t), I(t), R(t)) of system (1.5) with nonnegative initial value remain nonnegative and ultimately bounded. Firstly, by the first equation of system (1.5) we have Secondly, we show that I(t) is nonnegative for t ≥ 0. In fact, if there exists t 1 > 0 such that I(t 1 ) = 0 and I(t) > 0 for t ∈ (0, t 1 ), then by the second equation of system (1.5) we have I(t 1 + 1) = 0. This implies that I(t) ≥ 0 for all t ≥ 0.
Summarizing this analysis, we know that any solution of system (1.5) with nonnegative initial values remains nonnegative for all t ≥ 0, that is, all solutions of system (1.3) with nonnegative initial value remain nonnegative for all t ≥ 0.
In the following, we will show that the solutions of system (1.3) are ultimately bounded. It is clear that equation (1.4) has a unique equilibrium N * = Λ d , which is globally asymptotically stable, that is, lim t→∞ N(t) = N * . Therefore It follows that the omega limit set of system (1.3) is contained in the bounded feasible region Obviously, this region is positively invariant with respect to system (1.3), which implies that system (1.3) is ultimately bounded.
Since the first two equations of system (1.3) are independent of the third one, it suffices to consider the first two equations. Thus we replace h(I) with m in system (1.3) and restrict our attention to the following reduced model: The purpose of this paper is to show the dynamics of discrete system (1.6), In particular, we mainly focus on the various bifurcations that might occur around the equilibria. The paper is organized as follows. The saddle-node bifurcation, flip bifurcation, Neimark-Sacker bifurcation, and 1 : 1 resonance are shown in Sect. 2. In Sect. 3, a series of numerical simulations show that there are bifurcation and chaos in the discrete epidemic model. Finally, we give some conclusions.

Bifurcation analysis
Before bifurcation analysis, we first discuss the existence of the equilibria of system (1.6) and the global stability of the disease-free equilibrium. Letting R 0 = Λβ d(d+γ ) > 0 and A = mβ d(d+γ ) > 0, we have the following result. 1 , Proof To discuss the existence of equilibria of system (1.6), we need discuss the following equations: (2.1) When I = 0, we have m = 0. By the first equation of (2.1) we obtain S = Λ d . This implies that system (1.6) always has the disease-free equilibrium E 0 = ( Λ d , 0). When I > 0, we have m > 0. The first equation of (2.1) implies that S = Λ βI+d . Taking S = Λ βI+d in the second equation of of (2.1), we obtain the following equation: , that is, the following results hold: β , which implies that system (1.6) has one endemic equilibrium E 1 = (S 1 , I 1 ), where 1 , respectively. This means that system (1.6) has two endemic equilibria E 2 = (S 2 , I 2 ) and E 3 = (S 3 , I 3 ), where S i = Λ d+βI i , i = 2, 3, and Theorem 2.2 If R 0 < 1, then the disease-free equilibrium E 0 is globally asymptotically stable.
Proof We discuss the global stability of E 0 by constructing a Lyapunov function. Define Obviously, F is the mapping derived by the second It follows from Theorem 4.22 in [30] that E 0 is globally asymptotically stable.
In the following sections, we focus on bifurcations that may occur around the equilibria points.

Saddle-node bifurcation and bifurcation with 1 : 1 resonance
The Jacobian matrix at the endemic equilibrium E 1 = (S 1 , I 1 ) is The corresponding characteristic equation of J E 1 is which implies that h 1 (λ) = 0 has eigenvalues λ 1 = 1 and The following theorem is the case where the endemic equilibrium E 1 is a saddlenode bifurcation point when R 0 = 1 + A + 2 √ A. The existence of the endemic equilibria E 2 and E 3 is exhibited in Theorem 2.1, and the stability of both E 2 and E 3 is proved in the next section, Therefore we have the following: then system (1.6) will undergo a saddlenode bifurcation at E 1 . Moreover, when R 0 > 1+A +2 √ A, then system (1.6) has two endemic equilibria E 2 and E 3 .
In addition, that is, in this case, h 1 (λ) = 0 has the eigenvalues λ 1 = λ 2 = 1. The following theorem is the case where the endemic equilibrium E 1 is a 1 : 1 resonance point when 3 Λγ 2 . Then system (1.6) becomes (2.4) Furthermore, system (2.4) can be rewritten as Let The transformation for system (2.5) yields that where Using Sect. 9.5.2 in [31], we have It follows from Sect. 9.5.2 in [31] that system (1.6) will undergo a 1 : 1 resonance at E 1 In fact, the phenomenon of resonance with a double eigenvalue of 1 is similar to the Bogdanov-Takens bifurcation of the continuous system with double eigenvalues 0.

Flip bifurcation and Neimark-Sacker bifurcation
In this section, we first discuss the stability of the endemic equilibrium E 2 . Obviously, the Jacobian matrix at the endemic equilibrium E 2 = (S 2 , I 2 ) is as follows: and the corresponding characteristic equation of J E 2 is A, the endemic equilibrium E 2 is unstable, that is, we have the following: A, then the endemic equilibrium E 2 is unstable.
In the following, we discuss the stability of the endemic equilibrium E 3 and the bifurcation that may occur around E 3 when the stability of E 3 changes. Let It is clear that R 1 > 0 and R 2 > 0.
Proof To study the stability of E 3 , we give the Jacobian matrix at the endemic equilibrium E 3 = (S 3 , I 3 ): and the corresponding characteristic equation of J E 3 is It is clear that Therefore we have h 3 (-1) > 0 and h 3 (0) < 1 if 1 + A + R 1 < R 0 < 1 + A + R 2 , which implies that the endemic equilibrium E 3 is locally asymptotically stable, whereas if R 0 < 1 + A + R 1 or R 0 > 1 + A + R 2 , then the endemic equilibrium E 3 is unstable.
Because of 0 < d < 1, we have |λ 2 | = 1. The following theorem considers the case where the endemic equilibrium E 3 is a flip bifurcation point when R 0 = 1 + A + R 2 .

Moreover, the periodic-2 solution bifurcated from the endemic equilibrium E 3 is stable.
Proof In the case where R 0 = 1 + A + R 2 the endemic equilibrium is E 3 = (S 3 , I 3 ), where 3 , where the parameter τ is a new and dependent variable. The system (1.6) becomes: Using the transformation From the center manifold theory of discrete system we know that there exists a local manifold of system (2.8) [32]. The local manifold has the expansion 3 .
After substituting the expansion into system (2.8) and using the invariance property of the local manifold, the straightforward and careful calculation gives From the third equation of system (2.8) we have that z(t) is always constant. Therefore the one-dimensional model induced by the center manifold is 4 .
Proof Let x 1 (t) = S(t) -S 3 and y 1 (t) = I(t) -I 3 . Then system (1.6) becomes (2.9) The eigenvalues of the linearized matrix are Let R 0 be the bifurcation parameter, and let R * 0 = 1 + A + R 1 . A direct calculation yields that .
By a straightforward and tedious calculation we obtain that < 0, Using the Neimark-Sacker bifurcation theorem in [31], we obtain that there exists a Neimark-Sacker bifurcation when R 0 passes through R * 0 , which competes the proof.

Numerical simulation
To provide some numerical evidence for the qualitative dynamic behavior of the system (1.6), we will display some bifurcation diagrams to illustrate the analytical results and explore a new dynamics behavior as the parameters change.
We first show the stability of a disease-free equilibrium E 0 when R 0 > 1. Letting Λ = 0.1, d = 0.006, r = 0.1, m = 0.00001, and β = 0.0064, we have R 0 = 1.0063 > 1 and 1 + A + 2 √ A = 1.0202. Accordingly, we have 1 < R 0 < 1 + A + 2 √ A. In this case, Fig. 1(a) shows that the solutions of system (1.6) with different initial values tend to E 0 as t tends to infinity. Keeping the other parameters the same as in Fig. 1(a) except for β = 0.01, we have R 0 = 1.0535 > 1 and 1 + A + 2 √ A = 1.0206, that is, R 0 > 1 + A + 2 √ A. In this case, Fig. 1(b) shows that the solutions of system (1.6) also tend to E 0 as t tends to infinity. The numerical simulations illustrate that the disease-free equilibrium E 0 is stable when R 0 > 1.
Secondly, we show some examples about the existence of the endemic equilibrium E i , i = 1, 2, 3, Let Λ = 200, d = 0.07, r = 0.05, m = 0.1. Taking β as the variable parameter, we see in Fig. 2 the existence of the endemic equilibrium E i , i = 1, 2, 3. Then we display the saddle-node bifurcation around endemic equilibrium E 1 . Leaving the other parameters as in Fig. 2 Figure 3 shows that the solutions of system (1.6) with different initial values tend to E 3 as t tends to infinity, which implies that E 2 is unstable and E 3 is asymptotically stable, that is, E 1 is the saddle-node. shows that system (1.6) undergoes a process from periodic doubling to chaos. As β increases, firstly, the positive equilibrium loses its stability, and a stable two-cycle period appears, then the two-cycle period looses its stability, and a stable four-cycle period appears, and, finally, the period doubling process continues to chaos (see Fig. 7).

Discussion and conclusion
The discrete epidemic model is more suitable to describe the spread of diseases since the epidemiological data are usually collected in discrete time units, such as daily, weekly, or monthly. In this paper, we have studied the dynamical behavior of the discrete SIR epidemic model (1.3) with bilinear incidence rate and constant recovery by analyzing the discrete model (1.6). We have obtained sufficient conditions for the existence of the diseasefree equilibrium and multiple positive equilibria, the saddle-node bifurcation, the flip bifurcation, the 1 : 1 resonance, and the Neimark-Sacker bifurcation. Theoretic analysis shows that the disease-free equilibrium is globally asymptotically stable when the basic reproduction number is less than unity, and the numerical simulations illustrate that it is still asymptotically stable when the number is greater than unity. Model (1.3) is suitable to describe rotavirus, hand, foot, and mouth disease, influenza, and so on. These infectious diseases have the following common characteristics: (1) Most cases are mild; after isolation, some need take a symptomatic treatment, and some can recover without any treatment. (2) Only the tiny fractions of the cases have serious complications that need take treatment in hospital. Therefore, like rotavirus, hand, foot, and mouth disease, and influenza, the recovery form of the infective individuals consists of two parts, the spontaneous recovery and the recovery due to treatment. In this paper, we use γ I to describe the spontaneous recovery of the infective individuals without any treatment, and use h(I) to describe the recovery of the infective individuals with treatment. Naturally, the recovery function γ I + h(I) that we use here is a good description of the recovery of these diseases.
The dynamics of the discrete SIR epidemic model with bilinear incidence rate and constant recovery is much more complicated. We have investigated the stability of the diseasefree and positive equilibria and the bifurcation of the model analytically or numerically. In fact, the dynamic behavior of (1.3) has interesting practical significance.
Firstly, R 0 = Λβ d(d+γ ) is just the reproduction number of system (1.3) in the absence of the removal rate of the infective individuals. So, it is not the classic threshold of the epidemic dynamical models, and it cannot determines whether the disease will die out or not. However, 1 < R 0 < 1 + A + 2 √ A is a significant interval, which characterizes the importance of treatment for the disease control.
Secondly, the dynamic behavior of the model is different under different conditions when R 0 > 1+A +2 √ A. When 1+A + R 1 < R 0 < 1+A + R 2 , the endemic equilibrium E 3 is locally asymptotically stable. When R 0 = 1 + A + R 1 and 2d dγ (d+γ ) , the system will undergo a Neimark-Sacker bifurcation. When R 0 = 1 + A + R 2 , the system will undergo a flip bifurcation. Namely, in the case where the disease persists in population, the disease may tend to a certain equilibrium, or it may exist in the form of periodic shocks.
At last, the theoretical analysis show that in the case R 0 > 1 the disease-free equilibrium and endemic equilibrium not only coexist, but are both stable. This implies that the disease may still be extinct even if there is a stable endemic equilibrium, that is, if the initial state of the disease is within the stable domain of the disease free equilibrium, then the disease will disappear, whereas if the initial state of the disease is within the stable domain of endemic equilibrium, then the disease will persist.
There are still many challenging problems on the dynamics of the discrete system (1.3). The numerical simulations demonstrate that the disease-free equilibrium is stable when R 0 > 1. However, how to prove mathematically strictly the stability of the disease-free equilibrium E 0 when R 0 > 1 is a meaningful work. In addition, we also have no effective method to prove the global stability of the disease-free equilibrium E 0 when 1 < R 0 < 1 + A + 2 √ A. We think they may be related to the switching system theory and expect that some analytical results can be obtained on those problems in the future.
A summary of these calculations shows that R 1 ≥ 2 √ A.