Stability and bifurcations in a discrete-time epidemic model with vaccination and vital dynamics

Background The spread of infectious diseases is so important that changes the demography of the population. Therefore, prevention and intervention measures are essential to control and eliminate the disease. Among the drug and non-drug interventions, vaccination is a powerful strategy to preserve the population from infection. Mathematical models are useful to study the behavior of an infection when it enters a population and to investigate under which conditions it will be wiped out or continued. Results A discrete-time SIS epidemic model is introduced that includes a vaccination program. Some basic properties of this model are obtained; such as the equilibria and the basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_0$$\end{document}R0. Then the stability of the equilibria is given in terms of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_0$$\end{document}R0, and the bifurcations of the model are studied. By applying the forward Euler method on the continuous version of the model, a discretized model is obtained and analyzed. Conclusion It is proven that the disease-free equilibrium and endemic equilibrium are stable if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_0<1$$\end{document}R0<1 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_0>1$$\end{document}R0>1, respectively. Also, the disease-free equilibrium is globally stable when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_0\le 1$$\end{document}R0≤1. The system has a transcritical bifurcation when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {R}_0=1$$\end{document}R0=1 and it might also have period-doubling bifurcation. The sufficient conditions for the stability of equilibria in the discretized model are established. The numerical discussions verify the theoretical results.

discrete-time ones by difference equations. Recently, discrete models have gotten more interest because epidemic data are collected in discrete time intervals and numerical schemes also use discretization for solving differential equations. Moreover, the discretetime epidemic model exhibits more complex dynamics. Allen [3] studied some discretetime SI, SIR, and SIS epidemic models and showed that the simple discrete-time SI and SIR epidemic models without births or deaths mimic the behavior of the continuoustime models, while the behavior in the discrete-time SI, SIR, and SIS models with recovery or births differ from their continuous analogs. Castillo-Chavez and Yakubu [4] have investigated a discrete SIS model with complex dynamics. Brauer et al. [5] introduced a discrete-epidemic framework and highlight, emphasizing the final size of the epidemic, the similarities between single-outbreak comparable classical continuous-time epidemic models, and the discrete-time models. Farnoosh and Parsamanesh [6] investigated the stability and bifurcation in a discrete SIS model with bilinear incidence. The vaccination program implements not only on susceptible individuals but also on newcomers to the population. Parsamanesh and Mehrshad [7] performed a similar investigation on an SIS model with a temporary vaccination program with standard incidence. A discrete SIRS epidemic model with vaccination and a general infection probability function is investigated by Xiang et al. [8]. The vaccination performs only on susceptible individuals but not on the recruitment of the population. The local and global stabilities of disease-free equilibrium were derived as well as the local stability of the endemic equilibrium.
The discrete models have usually constructed either directly from properties of the disease and the population in which the propagation occurs, or by discretizing continuous model by a method based on numerical analysis such as the forward/backward Euler method, Mickens' non-standard method, or a mixed type formula which uses implicit and explicit Euler methods. For example, Roeger and Barnard [9] applied the central difference method to a continuous SIR epidemic model and showed the local stability of the equilibria. Liu et al. [10] presented four discrete epidemic models with the nonlinear incidence rate, using the forward Euler and backward Euler methods. They discussed the effect of two discretizations on the stability of the endemic equilibrium for these models. They concluded the forward Euler method makes the models have much richer dynamical behavior than the continuous models while the backward Euler method preserves the global asymptotic stability of endemic equilibria. Aranda et al. [11] provided a discrete model by discretization of the continuous-time model for transmission of Babesiosis disease in bovine and tick populations. They proved local stability of the equilibria and global stability of the disease-free equilibrium and obtained similar conclusions as they got in the continuous case. In particular, Mickens [12] consider the non-standard discretization for numerical methods which ensures positivity of the solutions of the difference equations. Recently, Izzo and Vecchio [13] obtained a discrete epidemic model by using a mixed type formula, which uses implicit and explicit Euler methods for discretization. Their discretization is similar to Mickens' non-standard one. But their discretization showed some good properties for positivity, boundedness, and global behavior of the solution. There are some works devoted to the local stability of a discrete model and few others investigated the global stability of the model. Hu et al. [14] established the criteria on the local stability of the disease-free equilibrium and endemic equilibrium for a class of SIRS epidemic model with a non-linear incidence, by using the linearization method and the expression of roots of the cubic polynomial equation. Ma et al. [15] studied the global stability of the endemic equilibrium of a discrete SIR model and get sufficient conditions by using the comparison principle. Cui and Zhang [16] gave sufficient conditions for the global stability of a discrete SIR model derived by employing a non-standard finite difference scheme. Van den Driessche and Yakubu [17] proved the global stability and persistence for some discrete epidemic models that are formulated for SEIR infections, cholera in humans, and anthrax in animals, under some demographic assumptions. Among these models, the susceptible-infected-susceptible (SIS) epidemic models are one of the well-known types of epidemic models. To consider the effect of vaccination as an efficient strategy to control and eliminate infections, it is possible to add a compartment for the vaccinated individuals to the SIS model and obtain the SIS epidemic model with vaccination, namely SIVS epidemic model [18,19]. These models may be deterministic [6] or stochastic [20], with constant [21] or variable [19] population size, and with standard [22] or bilinear [6,20] incidence rate. Motivated by the above studies, in this paper, we consider the discrete-time SIS epidemic model presented in [7], but with a different recruitment rate and without disease-caused death. Indeed, a model with the standard incidence in which a (perfect but temporary) vaccination program has been included. For this model using a similar approach used in [7], we investigate the local stability and bifurcations theoretically and numerically. Moreover, the global stability of the disease-free equilibrium of the model is proven. Furthermore, a generalization of the model is given by forwarding Euler discretization, and then its stability is studied. The effect of the vaccination in the model for controlling and eliminating the disease will also be shown.
The organization of the paper reads as follows: In the next section, the model is introduced, and equilibria of the model and its basic reproduction number are obtained. Two next sub-sections are devoted to studying the stability of the equilibria and bifurcations of the model, respectively. Then by using the forward Euler method, a discrete-time model is obtained from a continuous version of the model, and the stability of its equilibria is analyzed. After a numerical discussion, finally we summarize the results.

The model
Suppose that the individuals in a population are partitioned into susceptible individuals, infected individuals, and vaccinated individuals. Also, consider t as the appropriate time increment such that the changes in the model may take place at times 0, �t, 2�t, 3�t, . . . . The number of total individuals at time t = n t , for some n, is denoted by N t and numbers of individuals in other compartments in the same time are as S t , I t , and V t .
All possible changes in the model and transmissions between its sub-populations together with their transmission rates have been shown in Fig. 1. Here, all parameters are assumed to be nonnegative, and N and µ are positive. Also, µ is the natural death rate, β is the contact rate, γ is the cure rate, ǫ is the rate of losing immunity, while q and p are the vaccination rate in newcomers and susceptible individuals, respectively. The model can be illustrated by the following system of difference equations: The susceptible individuals become infected at standard incidence rate βS t I t /N t . Moreover, summing equations in system (1), we see that N t+1 = N t , and then the population size will remain a constant value. Thus, by letting V t = N − S t − I t , the corresponding difference equation is deleted and the following system of two difference equations is obtained: System (2) is considered under the following conditions, which are sufficient but not necessary for the nonnegativity of solutions.
These conditions are extracted from system (2) and they are natural requirements for the model. The first, state the rate at which the susceptible individuals who die or get infected or become vaccinated is less than one within a unit time. The second, says the rate at which the infected people who die or get recovered is less than one within a unit time.
The equilibria of the model are solutions of the following system: From the first equation, we must have either Ī = 0 or βS/N − (µ + γ ) = 0 . When Ī = 0 , the equilibrium is named the disease-free equilibrium and is written as (2) This equilibrium in which Ī � = 0 , is called the endemic equilibrium and is written as The quantity R 0 is referred to as the basic reproduction number of model (2) and is interpreted as the number of individuals who become infected by entering one infected individual into a fully susceptible population; see [23]. Thus it is reasonable intuitively, as we showed mathematically, that when R 0 < 1 , each infected individual transmits the infection to less than one other individual and it is expected that the infection vanishes. While in the case R 0 > 1 , each infected individual transmits the infection to more than one other individual and so the infection does not wipe out. Here, the assumption of R 0 was made directly from the expression of I * and condition of its positivity in endemic equilibrium. Also, there are some methods to calculate R 0 for discrete models; for example, see [24]. We see that R 0 is independent of the total population size N. Also and Therefore, we can state the following lemma about the existence of equilibria of the model. (2), the disease-free equilibrium Q 0 always exists and the endemic equilibrium Q * also exists if R 0 > 1.

Stability of the equilibria
We study stability of the system at an equilibrium by considering eigenvalues of the corresponding Jacobian matrix at that equilibrium. When eigenvalues are less than one, the system is stable.

Theorem 2
The disease-free equilibrium is stable if and only if R 0 < 1.

Proof
The Jacobian matrix of model (2) at (I, S) is Therefore, the Jacobian matrix at Q 0 is given by obviously, | 2 | < 1 by assumptions (3) and | 1 | < 1 if and only if R 0 < 1 .

Theorem 3
When R 0 > 1 , the endemic equilibrium Q * is stable and otherwise is unstable.
since µ + p + ǫ < 1 and βI * /N < 1 . Therefore, when R 0 > 1 , the Jury conditions are satisfied and the proof is completed. Now, in the next theorem we prove the global stability of the disease-free equilibrium.

Theorem 4
The disease-free equilibrium of model (2) is globally asymptotically stable if R 0 ≤ 1.
Proof First we divide the equations in system (2) by total population size N and letting s t = S t /N and i t = I t /N , we obtain the following equivalent system The basic reproduction number of this system is R eqv 0 = β µ+γ . Now, considering the function W (i t ) = i t , we see that W > 0 , and W = 0 if and only if i = 0 , that is at diseasefree equilibrium of (6). On the other hand, we have and thus when R eqv 0 ≤ 1 , we have W (i t+1 ) − W (i t ) ≤ 0 . Therefore W is a Lyapunov function [26] and the disease-free equilibrium of system (6) and as a result Q 0 is globally asymptotically stable.

Bifurcations of the model
In a discrete-time system, bifurcations occur at the equilibria of the given system when there exist some eigenvalues of the Jacobian matrix with module one. Indeed, for an eigenvalue , if = 1 , then a transcritical bifurcation (or flip bifurcation) occurs and when = −1 , a period-doubling bifurcation occurs; see [2,25]. While a Neimark-Sacker bifurcation, which is the same as the Hopf bifurcation in continuous systems [27], occurs if there is a pair of conjugate complex eigenvalues with module one, | | = 1.
As we saw, the eigenvalues of J (Q 0 ) are 1 = 1 − (µ + γ ) + (µ + γ )R 0 and 2 = 1 − (µ + p + ǫ) . Also, 1 = 1 if and only if R 0 = 1 and thus a transcritical bifurcation occurs at Q 0 when R 0 = 1 . On the other hand, 1 = −1 if and only if R 0 = 1 − 2 µ+γ , but this is impossible because µ + γ < 1 and R 0 becomes a negative value. This shows that a period-doubling bifurcation does not occur at Q 0 . Also, the eigenvalues of J (Q 0 ) are both real, and therefore a Neimark-Sacker bifurcation does not take place, too. Thus we have the following theorem.  Now, notice that as we concluded previously, P(−1) > 0 when R 0 > 1 . Also, R 0 = 1 implies βI * /N = 0 and this results in 2 − (µ + p + ǫ) = 0 which is impossible. These discussions state that a period-doubling bifurcation does not happen at Q * .
If we write the characteristic equation of J * as P( ) = 2 + a 1 + a 2 , we see that Hence, the roots of P( ) are both real and thus a Neimark-Sacker bifurcation cannot appear at Q * .
Remark 1 If we omit the restriction β < 1 from the system and allow β to take values greater than or equal to one, then from (7), we get Thus a period-doubling bifurcation occurs at Q * for β ≥ 1 if

The model obtained by the forward Euler discretization
The model described in Fig. 1 can be stated as a continuous-time model by the following system of ordinary differential equations (see [28,29]): We see that Ṅ = dN /dt = 0 and therefore the population size is constant. Similar to the discrete-time model, we get the following two-dimensional system by substituting V = N − S − I and omitting variable V from the system: Now in this section, we discretize and analyze model (9) by using the forward Euler method. Substituting Ṡ = (S t+1 − S t )/� and İ = (I t+1 − I t )/� , where is the fixed step size of the discretization, we obtain the discrete version of model as follows: It can be seen that the equilibria of this model and the corresponding basic reproduction number are similar to model (2). The disease-free equilibrium of discretized model, Q 0 d , always exists while, its endemic equilibrium, Q * d , exists only when R 0 > 1. .

Proof
The Jacobian matrix at endemic equilibrium is According to Jury conditions (Schur-Cohn criterion), the matrix J * d is stable (i.e., the roots of its characteristic equation P d ( ) = 2 + a 1 + a 2 lie inside the unit disk) if and only if the following conditions hold (see [30]): Here, a 1 = −tr(J * d ) and . Now, if we denote two roots as r 1 and r 2 (suppose r 1 < r 2 ), then P(−1) is positive when � < r 1 or � > r 2 , since b 2 > 0 . Moreover, we can easily see that b 1 − b 2 1 − 4b 2 > 0 , and thus r 1 > 0 . Therefore we can state, conditions (i)-(iii) hold if � < r 1 , because we also have r Remark 2 Model (2) was formulated straightly by considering a population and its transmissions. Also, the model can be concluded from discretized model (10) for = 1 and assumptions (3).

Discussion
In this section, we consider numerically theoretical results obtained in the paper. For this purpose assume that the parameters in the model are as q = 0.4 , p = 0.2 , γ = 0.15 , µ = 0.1 , and ǫ = 0.25 . Moreover, consider units of time and population as one day and one million individuals, respectively. Let the number of initial individuals in each subpopulations be as I 0 = 0.4 , S 0 = 0.8 , and V 0 = 0.5 . We take the contact rate β as the bifurcation parameter and get the bifurcation diagram as it is shown in Fig. 2. We see that at β = 0.443 , the dynamic of the system changes: The disease-free equilibrium that was stable for values β < 0.443 becomes unstable and instead the endemic equilibrium becomes stable. Indeed, at β = 0.4435 , we have R 0 = 1 and a transcritical bifurcation occurs. Moreover, it is seen that at β = 2.428 , the endemic equilibrium becomes unstable and a period-doubling bifurcation happens and after that, the system remains unstable. This value for β also is obtained according to Remark 1 as β = 2.4279 . Figure 3 shows the Lyapunov exponents of the Jacobian matrix for the same values of β . Here also it is observable that for values β = 0.443, 2.428 , and 3.235, the Lyapunov exponent is positive as seen in the bifurcation diagram. Figure 4 presents solutions of the system for various values of β and the behavior of the solutions is the same as we expect from the bifurcation diagram and the Lyapunov exponents. For β = 0.4 , we have R 0 = 0.9018 < 1 and as we expect from Theorem 2, the disease will vanish. While, for values β = 0.55 , β = 2.45 , and β = 2.5 , we have R 0 = 1.2400 , R 0 = 5.5236 , and R 0 = 5.6364 , respectively, that all  Bifurcation diagram for I t in terms of in model (10) are greater than one and according to Theorem 3, the infection remains at a positive level. In addition, Fig. 5 displays some parts of solutions of infected population I t for different values of β . It is observable that the behavior of solutions corresponds to those are in the bifurcation diagram. The effect of the discretization of the continuous-time model (9) by applying the forward Euler method has been considered in Fig. 6. The bifurcation diagram shows the dynamics of the infected population in the discretized model (10) when the step size varies. The contact rate has been supposed as β = 2.7 and other parameters are the same as preceding simulations. As it was established in Theorem 8, the endemic equilibrium Q * d is stable for � < � * = 0.8946 while for greater values of , the endemic equilibrium becomes unstable.

Conclusions
In this paper, we introduced and studied an SIS epidemic model that includes a vaccination program. The equilibria of the model were detected: The disease-free equilibrium Q 0 in which the infection will be extinct, and the endemic equilibrium Q * in which the disease will persist in the population. It was proved that under some assumptions on parameters for positivity of solutions, Q 0 and Q * are stable if R 0 < 1 and R 0 > 1 , respectively. It was proven Q 0 is also globally asymptotically stable when R 0 ≤ 1 . Thus the basic reproduction number R 0 plays important role in determining the dynamics of the model. To clear the effect of vaccination, a sensitivity analysis is performed on R 0 to determine how sensitive the model is to changes in the value of parameters that are related to vaccination. The normalized forward sensitivity index of R 0 concerning parameter x has been defined in [31] as � R 0 x = x R 0 × ∂R 0 ∂x . The sensitivity indices for vaccination proportions p and q are calculated as � R 0 q = −qµ (1−q)µ+ε < 0 and � R 0 p = −p µ+p+ε < 0 . These indices are both negative meanings R 0 is inversely related with p and q; an increase in parameters will cause a decrease in value of R 0 , that is more vaccinated individuals, the infection becomes more controllable. On the other hands, The basic reproduction number for the corresponding model without vaccination (i.e. p = q = 0 ) is R 0 = β µ+γ . Thus R 0 = (1−q)µ+ε µ+p+ε R 0 and R 0 <R 0 . When R 0 > 1 , then the vaccination must be used such that R 0 < 1 to eliminate an infection. Furthermore, the bifurcations of the model were investigated and it was proved that when R 0 = 1 system has a transcritical bifurcation and although the Neimark-Sacker bifurcation does not appear, it may have a period-doubling bifurcation if we ignore the restriction β < 1 . To study the discretization of the continuous version of the model, we applied the forward Euler method and analyzed the effect of step size of the discretization on the dynamics of the model. We established the sufficient condition for stability of disease-free equilibrium Q 0 d and endemic equilibrium Q * d in the discretized model. Finally, we examine the results obtained in the paper in numerical example by considering the bifurcation diagram, the Lyapunov exponents of the Jacobian matrix, and graphs of solutions for values of β and . It was observed that the numerical discussions verify the theoretical results.