Global stability of the endemic equilibrium of a discrete SIR epidemic model

The basic reproductive number R0 of a discrete SIR epidemic model is defined and the dynamical behavior of the model is studied. It is proved that the disease free equilibrium is globally asymptotically stable if R0<1, and the persistence of the model is obtained when R0>1. The main attention is paid to the global stability of the endemic equilibrium. Sufficient conditions for the global stability of the endemic equilibrium are established by using the comparison principle. Numerical simulations are done to show our theoretical results and to demonstrate the complicated dynamics of the model. Electronic supplementary material The online version of this article (doi:10.1186/1687-1847-2013-42) contains supplementary material, which is available to authorized users.


Introduction
Mathematical models have played a significant role in describing the dynamical evolution of infectious diseases. SIR model is one of the classical epidemic models with compartment structure, suitable for the transmission of infectious diseases that confer longlasting immunity such as chicken pox and SARS. The host population is divided into three epidemiological groups: the susceptibles, the infectives, and the removed/recovered. The transmission dynamics of an infectious disease is described by modeling the population movements among those epidemiological compartments.
There is an increasing interest in the study and application of discrete epidemic models. Allen et al. have studied some discrete-time SI, SIR, and SIS epidemic models, compared the dynamics of deterministic and stochastic discrete-time epidemic models, and also applied the discrete epidemic model to the spread of rabies [-]. Allen and Driessche have given the basic reproductive number for some discrete-time epidemic models []. Castillo-Chavez and Yakubu have investigated a discrete SIS model with complex dynamics []. Zhou et al. have formulated and discussed the dynamical behavior of agestructured SIS models [, ]. A stage-structured model, a discrete epidemic model with seasonal variation in environment, a discrete tuberculosis transmission model, and many other discrete epidemic models have been constructed, studied, and applied [-]. One reason for the upsurge of discrete epidemic models is that discrete models have advantages in describing an infectious disease since epidemic data are usually collected in discrete time units, which would make it more convenient to use discrete-time models [].
The studies on discrete epidemic models are relatively few compared with those on continuous epidemic models. Most research works on discrete epidemic models concern the http://www.advancesindifferenceequations.com/content/2013/1/42 definition of the basic reproductive number, the global stability of the disease free equilibrium, the persistence of diseases, the existence and local stability of endemic equilibria, the existence of flip bifurcation and Hopf bifurcation. The results on the global stability of the endemic equilibrium are quite few for discrete epidemic models due to the following two factors: () There is not enough effective theory and methods for global stability of discrete dynamical systems; () The discrete model exhibits much more complicated dynamical behaviors than its corresponding continuous model.
We will study the global stability of the endemic equilibrium of a discrete SIR model and get sufficient conditions. The discrete SIR model, the biological requirement of the model, the basic reproductive number, and the invariant domain of the model are given in the next section. The global stability of disease free equilibrium and the persistence are discussed in Section . The sufficient conditions for the global stability of endemic equilibrium are obtained in Section . Numerical simulations and application are presented to demonstrate our theoretical results, to show the complex dynamics and application example of the model in Section . Conclusions and discussions are included in the last section.

The discrete SIR model
There are different approaches to model infectious disease evolution in discrete time. The recurrent difference equations from the discretization of continuous differential equation models is one of the direct and frequently-used modeling approaches. This kind of model can be well understood in application under reasonable assumptions though there are some limitations on the range of parameters. We will adopt this approach to formulate model and assume that the infected people will obtain permanent immunity after they get recovered from infection.
Let N(t) denote the number of the host population at time t. According to the disease transmission mechanism, the host population is grouped into three epidemiological compartments. Let S(t), I(t), and R(t) be the number of individuals in the susceptible, infective, and removed/recovered compartments at time t, respectively. In addition to the death and recruitment, there are population movements among those three epidemiological compartments from time unit t to time t + . We assume that the recruited individuals (by birth and immigration) are constant and enter the susceptible compartment. After one unit time, the susceptible individuals may stay in the susceptible compartment, or get infected and move to the infectious compartment, or die. The individuals in the infective compartment can keep being the infective, or get recover and transfer to the recovered compartment, or die. The individuals in the recovered compartment never leave the compartment unless they die. The discrete SIR model is where is the constant recruitment rate of the population, β is the transmission rate, μ is the natural death rate, δ is the death rate caused by disease, and γ is the recovery rate. The biological background requires that all parameters be non-negative.
There is abundant amount of research into SIR models since they capture the basic evolution mechanism of the infectious diseases when the recovered individuals will acquire http://www.advancesindifferenceequations.com/content/2013/1/42 life-long immunity. A lot of results on the existence and global stability of the endemic equilibrium have been obtained for continuous SIR models with various transmission rate. However, the result on the global stability of the endemic equilibrium for discrete SIR model is rare, even for the model with mass action incidence. We will try to give sufficient conditions for the global stability of the endemic equilibrium. Those results on global stability of the endemic equilibrium can promote the study on this challengeable problem.
Adding all equations in (), we see that the number of host population, N(t), satisfies The equationÑ(t + ) = + (μ)Ñ(t) has a unique equilibriumÑ * = μ , which is globally asymptotically stable, namely, lim t→∞Ñ (t) =Ñ * . The comparison principle implies that In the analysis of model (), we assume that The epidemiological interpretation requires the solution of model () with initial values satisfying () to be non-negative. This requirement can be met if the following two inequalities hold: Those two conditions in () are natural requirements for model (). The first inequality, β + μ < , says that the percentage of the susceptible individuals who die or get infected is less than one within a unit time. The second inequality, μ + δ + γ < , states that the percentage of the infected people who die or get recovered is less than one within a unit time. It is easy to verify that those two inequalities ensure S(t) ≥ , I(t) ≥ , and R(t) ≥  for all t ≥  if the initial values satisfy ().
In the rest of the paper, we assume that the initial conditions and the parameters satisfy () and (), respectively. It is easy to verify that the domain is a compact, positively invariant set for model ().  is also a global compact attractor of system () since it attracts all positive orbits with an initial value (S  , I  , R  ) ∈ R  + . The basic reproductive number, R  , is fundamental and widely used in the study of epidemiological models. R  is the average number of secondary infections produced by one infected individual during the entire course of infection in a completely susceptible population. R  often serves as a threshold parameter that predicts whether an infection dies out or keeps persistence in a population. Following the idea and framework given in [, , ], we obtain the formula of the basic reproduction number of model (), R  = β μ+δ+γ . The magnitude of R  plays a crucial role in determining the dynamical behavior of model (). http://www.advancesindifferenceequations.com/content/2013/1/42

Extinction and persistence for the disease
This section focuses on the disease extinction or persistence, which is determined by the stability of the disease free equilibrium and the existence of endemic equilibrium of model (). It is obvious that E  ( μ , , ) is an equilibrium of model (). E  is called the disease free equilibrium since I and R components are zero. The stability of the disease free equilibrium E  is given in the following theorem.
Proof The linearized matrix of () at the disease free equilibrium E  is Obviously, the eigenvalues of matrix A are λ  = λ  = μ and λ  =  + βμδγ , respectively. The conditions β + μ < , μ + δ + γ < , and R  <  imply that |λ i | < , i = , , . Therefore, the disease free equilibrium E  is asymptotically stable. When R  > , we have From the second equation of model () and the fact S(t) ≤ N(t), we have The conditions μ + δ + γ <  and R  <  imply that  <  -(μ + δ + γ )( -R  ) < . The recurrent use of inequality () yields The inequality () implies that lim t→∞ I(t) = . From the fact lim t→∞ I(t) = , we know that for any ε > , there exists a large positive integer T  such that I(t) < ε when t ≥ T  . Consequently, The equationR(t +) =R(t)+γ ε -μR(t) has a unique equilibriumR * = γ ε μ , which is globally asymptotically stable. The comparison principle implies that there exists an integer T  > T  such that R(t) ≤R(t) < γ ε μ if t > T  . The arbitrary ε implies that lim t→∞ R(t) = . From () we have From the left inequality of () and the comparison principle, we know that for any given ε  > , there exists an integer T  > T  such that N(t) ≥ -δε με  for all t > T  . From the right inequality of () and the comparison principle, we know that for any given ε  > , http://www.advancesindifferenceequations.com/content/2013/1/42 and the arbitrary ε, ε  , and ε  imply that lim t→∞ N(t) = μ , i.e., lim t→∞ S(t) = lim t→∞ (N(t)- imply that the disease free equilibrium of () is globally asymptotically stable when R  < .
R  >  means that the average number of a new infection by an infected individual is more than one. The epidemiological interpretation indicates that the disease may keep persistent in the population. The following theorem confirms the persistence of the disease in the case R  > .
Theorem  If R  > , then the disease will keep persistence in the population, i.e., there exists a positive ε such that the solution of model has a positive equilibrium S * = μ , which is globally attractive. By using Lemma . in has a unique equilibrium N *  = μ+δ , and N(t + ) = N(t) + -μN(t) has a unique equilibrium N *  = μ , which is globally asymptotically stable. Therefore, for any given ε > , there exists a T  >  such that (μ+δ)ε ≤ N(t) ≤ μ + ε for all t ≥ T  .
If R  > , then we can prove that there exists a small positive number σ such that If the conclusion in () does not hold, then for any positive number σ , there exist a point The inequality in () implies that When t > T  , the equations in () imply that () From the first inequality in (), we know that there exists a number T  > T  such that N(t) ≤ μ holds for all t > T  . When t > T  , we substitute the inequality N(t) ≤ μ into the second inequality of () to obtain If we choose σ small enough, then the condition R  >  implies that From the inequalities in () and (), we have that lim t→∞ I(t) = ∞. The limit lim t→∞ I(t) = ∞ contradicts the inequality I(t) < σ in (). The contradiction comes from the assumption given in (). Therefore, the conclusion in () holds true. Subsequently, W s (  ) ∩ X  = ∅ and  is isolated in X. From Theorem .. and Remark .. in [], it follows that is uniformly persistent with respect to (X  , ∂X  ). Furthermore, Theorem .. in [] implies that the solutions of model () are uniformly persistent with respect to (X  , The conditions in Theorem  and Theorem  are very simple though the proof is long and complicated. It is very easy to determine whether the disease goes to extinction or keeps persistence in the population by the magnitude of the basic reproductive number, R  . The underlining mechanism of the nice threshold result given in Theorem  is that the number of the infected population will gradually become lower and lower if the new infection created by an individual during his/her infection period is less than one. The epidemiological interpretation of Theorem  says that there is at least a certain number of infected population if the new infection created by an individual during his/her infection period is more than one. The necessary condition is to eliminate an infectious disease, to reduce R  , and to have R  < . The basis reproductive number, R  , is a very useful quantity in mathematical analysis and epidemiological application. http://www.advancesindifferenceequations.com/content/2013/1/42

The stability of the endemic equilibrium
This section deals with the existence and stability of the endemic equilibrium of model (). Let E * (S * , I * , R * ) be the endemic equilibrium of model (). Then S * , I * , and R * satisfy the following equations: where N * = S * + I * + R * . By straightforward and careful calculations, we know that model () has a unique endemic equilibrium when R  >  with , The local stability of the endemic equilibrium E * is given in the following theorem.
Proof In order to discuss the stability of the endemic equilibrium of model (), we consider the following equivalent system: When R  > , the positive equilibrium of system () is (N * , I * , R * ), where .
The linearization matrix of () at the positive equilibrium (N * , I * , R * ) is The characteristic equation of matrix A  is We see that the equation φ(λ) =  has an eigenvalue  < λ  = μ < . Therefore, in order to determine the stability of the positive equilibrium of model (), we discuss the roots http://www.advancesindifferenceequations.com/content/2013/1/42 of the following equation: When R  > , the calculation yields Furthermore, the constant term satisfies The Jury criterion [] implies that the two roots, λ  and λ  , of the equation φ  (λ) =  satisfy |λ  | <  and |λ  | < . The linearization theory implies that the positive equilibrium (N * , I * , R * ) of system () is locally asymptotically stable if R  > , i.e., the endemic equilibrium E * of system () is locally asymptotically stable.
The global stability of the endemic equilibrium E * of model () is quite difficult. Sufficient conditions are presented in two theorems for the special case δ =  and for the general case.
Proof When δ = , the host population N(t) = S(t) + I(t) + R(t) satisfies which has the solution From the solution of the host population, we have that lim t→∞ N(t) = N * = μ , and the limit of N(t) is independent of the initial value N  . The global stability of the endemic equilibrium, E * , of model () is equivalent to the global stability of the endemic equilibrium E * N (N * , S * , I * ) of the following model: () The first equation of model () is independent of other two state variables S(t) and I(t).
The fact that lim t→∞ N(t) = N * leads to the following limit model: From those two inequalities in () and the comparison theorem [], we know that for any small ε > , there exists a positive integer T l such that L l

I(t + ) ≥ I(t) + βμ(L l  -I(t))I(t) -(μ + γ )I(t).
() We consider the following two difference equations corresponding to the inequalities in (): () http://www.advancesindifferenceequations.com/content/2013/1/42 We substitute I m  (t) = +βμL m  -(μ+γ ) βμ x m  (t) and I l  (t) = +βμL l  -(μ+γ ) βμ x l  (t) into the first and the second equations of () to have The variables x m  (t) and x l  (t) satisfy the quadratic difference equation. The well-known result of the famous population model x(t + ) = rx(t)(x(t)) says that x *  =  is the unique and globally asymptotically stable equilibrium if  < r < , whereas x *  =  - r is the unique positive equilibrium if  < r <  and x *  is globally asymptotically stable. The result on the quadratic population model x(t + ) = rx(t)(x(t)) implies that the first equation in () has a positive equilibrium I m  * = L m  -μR  , which is globally asymptotically stable. A similar argument implies that the second equation in () has a positive equilibrium I l  * = L l  -μR  , which is globally asymptotically stable if R  > μ+γ μ . The asymptotical stability and the comparison theory imply that there exists a T i ≥ T l such that I l  ≤ I(t) ≤ I m  for all t > T i , where I l  = I l  *ε and I m  = I m  * + ε. When t ≥ T i , we substitute the inequality I l  ≤ I(t) ≤ I m  into the first equation of () and have () From () and a similar argument, we can obtain that there exists a positive integer T l When t ≥ T l , the inequality L l  ≤ L(t) ≤ L m  and the second equation of () imply that

I(t + ) ≥ I(t) + βμ(L l  -I(t))I(t) -(μ + γ )I(t). ()
A similar argument implies that there exists a T i ≥ T l such that I l Equations in () hold when t > T i and I l  ≤ I(t) ≤ I m  . From the mathematical induction, we know that there exist sequences T kl , T ki , L l k , L m k , I l k , and I m k such that I l k < I(t) < I m k for all t > T ki . Furthermore, I l k and I m k satisfy the following recurrence equations: () http://www.advancesindifferenceequations.com/content/2013/1/42 () is a linear system of difference equations. System () has a positive equilibrium P i * (I l * (ε), I m * (ε)), where Let λ  and λ  be two roots of the linearized matrix of system () at the equilibrium P i * . It is easy to see that The stability theory of difference equations implies that the equilibrium P i * of () is globally asymptotically stable, i.e., lim k→∞ I l k = I l * (ε) and lim k→∞ I m k = I m * (ε). From the expressions of I l * (ε) and I m * (ε), we have that lim ε→ I l . From the inequality I l k < I(t) < I m k and those limits, we obtain that lim t→∞ Similarly, we can prove that the sequences {L l k } and {L m k } satisfy a linear system of difference equations. The sequences {L l k } and {L m k } tend to the same limit when k tends to infinity and ε tends to zero, i.e., The inequality L l k < L(t) < L m k and the limit imply lim t→∞ L(t) = (βμ+γ (μ+γ )) Therefore, the endemic equilibrium of system () is globally asymptotic stable when R  > μ+γ μ and γ < μ. Equivalently, the endemic equilibrium of system () is globally asymptotic stable.
The comparison principle is the main idea to prove Theorem . The limitation of the method and the construction of the comparison equations may lead to the imposed conditions R  > μ+γ μ and γ ≤ μ. We do hope that the global stability of the endemic equilibrium can be proved under the condition R  > . When δ > , the global stability condition of endemic equilibrium is more complicated. We use the same idea to get the sufficient stability conditions.

Proof Let L(t) = S(t) + I(t) and N(t) = L(t) + R(t).
We consider the equivalent model

N(t + ) = N(t) + -μN(t) -δI(t),
() From () and the comparison theorem, we know that for any small ε > , there exists a positive integer T n such that N l From the second equation of system (), we have The comparison equations corresponding to those inequalities in (), When t ≥ T i , we use the inequality I l  ≤ I(t) ≤ I m  in the first two equations of system () to get A similar argument implies that there exists a positive integer T c such that N l When t > T c , we use those estimates of N(t) and L(t) in the third equation of () and obtain the following inequalities: () When μ  > (μ+δ)(μ+δ+γ )  +μ+δ+γ and μ+δ+γ μ < R  < ( +  μ+δ+γ ) μ μ+δ , a similar procedure as aforementioned can imply that there exists an integer T i > T c such that I l By using the mathematical induction, we obtain the sequences T ki , N l k , N m k , L l k , L m k , I l k , and I m k such that Furthermore, I l k and I m k satisfy the following recurrence equations: If the conditions of Theorem  hold, then we know that the two eigenvalues of matrix B satisfy |λ j | <  (j = , ), and the matrix series ∞ k= B k converges to (I -B) - . Under the http://www.advancesindifferenceequations.com/content/2013/1/42 conditions of Theorem , we can have lim k→∞ z(k) = (I -B) - b, i.e., z * = (I m * , I l * ) = ((I -B) - b) τ is the globally stable equilibrium of (). Further calculation shows that After taking ε → , we have and A similar argument implies that .
Those limits lead to Therefore, the endemic equilibrium of system () is globally asymptotic stable when the conditions of Theorem  hold.
There are some parameter values which can satisfy the conditions of Theorem . For example, if δ = γ = , those conditions become μ > ,  + μ > ,  < R  <  +  μ . Then, for μ > , we know that the conditions of Theorem  will hold if δ and γ are small enough.

Numerical simulations
In this section, we carry out numerical simulations to demonstrate our theoretical results and the complex dynamics of model () for several sets of parameters and initial values. We use the SIR model to simulate the mumps infection in China, too.

Simulations of model (1)
The simulations in this subsection demonstrate our theoretical results or the complex dynamics of the model. The parameter values for those simulations are selected mathematically to let the model exhibit required dynamics.
When δ = , we choose μ = ., γ = ., = , and β = ., then γ < μ, μ+γ μ = ., and R  = . > .. Theorem  implies that the endemic equilibrium E * (S * , I * , R * ) http://www.advancesindifferenceequations.com/content/2013/1/42 of () is globally asymptotically stable. The number of the infectious individuals for two solutions with different initial values are shown in Figure (a). From Figure (a) we see that the number of infected people approaches its equilibrium value as t tends to infinity. For the same values of δ, μ, γ , and , the straightforward calculation yields R  ∈ (, .) if β ∈ (., .), which does not satisfy the condition R  > μ+γ μ of Theorem . On the other hand, the simulation shows that the endemic equilibrium of model () may be globally asymptotically stable even if the condition R  > μ+γ μ of Theorem  does not hold. The solution curves of the infected people for β = . and two different initial values are given in Figure (b). Those simulations remind us that the condition R  > μ+γ μ of Theorem  may not be necessary for the global stability. It is a pity that we cannot prove the global stability without the condition R  > μ+γ μ . For the case δ > , we choose δ = ., μ = ., γ = ., and = . Then we have max{ δ μ-δ-γ , μ+δ+γ μ } = ., ( +  μ+δ+γ ) μ μ+δ = ., μ  > (μ+δ)(μ+δ+γ )  +μ+δ+γ , and μ > δ + γ . Theorem  implies that the endemic equilibrium E * is globally asymptotically stable when . < R  < .. Figure   From the epidemiological interpretation, we assume that the inequalities  < μ+δ +γ <  and β + μ <  hold. If those two inequalities hold, then all solutions of model () with positive initial values are non-negative, and the numerical simulations show that model () does not exhibit complex dynamics. If those two inequalities do not hold, then model () may exhibit much more complicated dynamical behaviors. The numerical simulation demonstrates that there exists a sequence of period-doubling bifurcation to chaos. Let us take μ = ., γ = ., δ = ., and =  to investigate the period doubling process of model (). For small R  (a linear function of β), the endemic equilibrium is unique and globally asymptotically stable (see Figure (a)). When R  passes through ., the endemic equilibrium losses its stability, and a stable periodic solution of period  appears (see Figure (b)). When R  passes through ., the periodic solution of period  losses its stability, and a stable periodic solution of period four appears (see Figure (c)). As R  increases further, the periodic solution of period  becomes unstable, and a periodic solution of period  appears. The numerical simulation shows that the period-doubling bifurcation may continue and go to chaos (see Figure (d)).

Application of model (1)
The simulation in this subsection is an application of the SIR model to the mumps infection in China. The model parameters are estimated according to the demographic and epidemiological data in China from  to . We have compared the model simulation result of the number of annual mumps cases to the reported number of notifiable diseases. http://www.advancesindifferenceequations.com/content/2013/1/42 Mumps, caused by the mumps virus, is an acute respirator infectious disease and it is spread from person to person by coughing or sneezing. Mumps usually starts with a fever and headache for a day or two. It then presents with swelling and soreness of the parotid salivary gland. The main symptom of mumps is the swollen, painful salivary glands on one or both sides of the face. The entire course of mumps infection is  to  days. The infected people will obtain permanent immunity after they get recovered from infection. The SIR model is suitable to describe the infection of mumps due to its permanent immunity. The host population is chosen to be individuals who are younger than  since mumps is common in children and adolescents, but rare among adults. The values of model parameters are estimated from the demographic and epidemiological data in China from  to .
In the application of the SIR model, the simulation time unit is taken as one day. All the parameter values will be estimated on the basis of that time unit. From the statistical yearbook of China, we know that the number of the host population (less than twenty years old) in  is  million []. The annual birth rate in China from  to  is approximately /,, and the recruitment rate is taken as = , [-]. The annual death rate of the population is ./, [, ], hence, the daily death rate is ./,,. Since the age of the host population is less than twenty, an individual should be removed when his/her age is over twenty. The daily removed rate is taken to be /,, and the value of parameter μ is estimated to be μ = . + . = .. From the fact that the entire course of mumps infection is  to  days, we take γ = / = . []. From the report on notifiable infectious diseases in China, we take δ = ..
The initial time is taken to be January first of . From the data of China population and notifiable diseases in  [], the initial values are taken to be S() = ,,, I() = ,, and R() = ,,, respectively. By using the method of parameter estimation in [], we have β = .. After the parameters and initial values are set, we apply the SIS model to simulate the mumps infection in China from  to . The simulation results are shown in Figure . The top left plot is the number of the susceptibles, S(t), which shows a slow downward tendency. From the demographic data in China, we know that the number of the population younger than twenty is  million in ,  million in , and  million in  [, , ]. The change of the host population in our simulation is consistent with the observed data. The top right is the number of the infective, I(t), which shows the slow upward tendency. The bottom left is the number of the recovered R(t). The curve in the bottom right plot is the number of new mumps infection cases from  to  by a model simulation. The piecewise curve with stars is the number of annually reported mumps cases out of notifiable infectious diseases from  to  []. The comparison between the model simulation and the notifiable disease report demonstrates that our model can give a good prediction for infection diseases.

Conclusion
SIR models are suitable to describe the transmission of infectious diseases with life long immunity. A lot of continuous SIR models with various transmission rates have been formulated and studied. The global stability of the endemic equilibrium of the continuous SIR model has been investigated extensively and many sufficient conditions have been obtained. Nevertheless, the discrete SIR model and its dynamics are quite few. In this paper, http://www.advancesindifferenceequations.com/content/2013/1/42 we have formulated a discrete SIR epidemic model and studied its asymptotical behaviors. Especially, we have got some sufficient conditions for the global stability of the endemic equilibrium. We have also applied the SIR model to the mumps infection in China. The model simulation results match the reported data of the notifiable diseases well.
Although the SIR model considered in this paper is quite simple, it exhibits very complicated dynamical behavior. The complicated dynamics of the simple SIR model reminds us that we cannot expect to have the global stability of the endemic equilibrium when R  > . Other conditions should be given for the global stability. In this paper we have obtained sufficient conditions for the global stability of the endemic equilibrium. Although those conditions given in our paper are not satisfactory, our result is a good exploration of the challengeable problem. We expect to prove the global stability of the endemic equilibrium when R  is not large. The comparison principle used in this paper and the Lyapunov function may help us to have better results.