A Nonstandard Finite Difference Scheme for SIS Epidemic Model with Delay : Stability and Bifurcation Analysis

A numerical scheme for a SIS epidemic model with a delay is constructed by applying a nonstandard finite difference (NSFD) method. The dynamics of the obtained discrete system is investigated. First we show that the discrete system has equilibria which are exactly the same as those of continuous model. By studying the distribution of the roots of the characteristics equations related to the linearized system, we can provide the stable regions in the appropriate parameter plane. It is shown that the conditions for those equilibria to be asymptotically stable are consistent with the continuous model for any size of numerical time-step. Furthermore, we also establish the existence of Neimark-Sacker bifurcation (also called Hopf bifurcation for map) which is controlled by the time delay. The analytical results are confirmed by some numerical simulations.


Introduction
In this paper we reconsider a SIS epidemic model with maturation delay developed in [1]: Here is a birth rate function.The size of mature population (2.c) Note that Equation (2.c) is decoupled and is called the Nicholson's blowflies equation which proposed by [2].The dynamics of Equation (2.c) have been studied by many authors; see e.g.[3,4].
In [1], the basic reproduction number for system (2) has been identified as It has also been shown in [1] that 0 determines the existence of equilibria as well as their stability properties.They have shown numerically that if the delay R  is sufficiently large then the positive solutions of system (2) oscillate about the positive equilibrium.The existence and stability of Hopf bifurcation were established by Wei and Zou [5].In this case, the bifurcation is controlled by parameter p.The bifurcation analysis of system (2) using time delay as the control parameter has been done by Chi, Qu and Wei [6].
For practical purposes, we need to do numerical simulations and therefore we have to transform the continuous model (2) into a discrete system.We expect that the dynamical properties of the discrete system are in accordance with its continuous counterpart (2).Kunnawuttipreechachan [7] and the author [8] considered the Euler discretization of system (2) and showed that the discrete system has exactly the same equilibria as those of system (2).Using different method of analysis, Kunnawuttipreechachan [7] and the author [8] derived the sufficient conditions of the numerical step-size for equilibria to be asymptotically stable.It was concluded that the stability conditions of equilibria of the discrete system obtained by the Euler method are consistent with system (2) only if the numerical time-step is relatively small.
To overcome the dependence of stability condition on the time-step size, we will apply a nonstandard finite difference (NSFD) scheme.This method, which is developed by Mickens [9,10], has been applied to various problems; see e.g.[11][12][13][14][15][16], in which the numerical solutions preserve dynamical properties of the continuous model.We will show that the discrete SIS epidemic model with a delay obtained by the NSFD method maintains the stability properties of equilibria irrespective of the size of numerical time-step.Besides the stability conditions, the existence of bifurcation of the discrete model will also be investigated.

Discrete SIS Epidemic Model with Delay
Using the fact that      , we will only consider the last two equations of system (2).Using a transformation   To discretize system (4) we first consider the step-size of the form 1 h k where k is a positive integer.Then, applying the forward difference scheme for the derivative and a nonlocal approximation for the right hand sides of system (4) yields a system of difference equations where j I and j N are numerical approximation of , respectively.The numerical scheme (5) is a nonstandard because it uses a nonlocal approximation; see Mickens [9,10].Here we consider initial conditions 0,1, 2, j     0 0 where . It is easy to show that the implicit scheme (5) can be arranged to get its explicit version, i.e.
Direct calculations show that the discrete system (7) has a unique equilibrium:   . This equilibrium is called the disease free equilibrium (DFE).On the other hand, if 0 and 1 R     then, in addition to the DFE, there also exists an endemic equilibrium (EE):  , We observe that those equilibria are exactly the same as those of the continuous system (4); see e.g.[7].

Stability and Neimark-Sacker Bifurcation Analysis
It is well known that the stability of equilibrium of a dynamical system depends on the distribution of the zeros of its associated characteristics equation.In this section, the distribution of the roots of the characteristics equation will be analyzed using the following results of Zhang, Zu and Zheng [17].
Theorem 1 (Zhang, Zu and Zheng [17]).Suppose that is a bounded, closed and connected set; and  is a parameter.Then as  varies, the sum of the order of the zeros of   can change only if a zero appears on or crosses the unit circle. (5)

Disease Free Equilibrium
First we perform a linearization of system (7) about the The linearized system around the DFE is given by By introducing new variables we can rewrite system (8) in the form where and the constant matrix A is defined by where .
  and   then there exists a   such that all roots of Equation (10) have modulus less than one for 0      .
Proof.The trivial root of Equation ( 10) is Other roots of Equation (10) are determined by whenever p   .Consequently all roots of Equation ( 10) lie in 1   for all sufficiently small 0   , and the existence of the maximal   follows.□ then Equation (11) has no root with modulus one for all 0   .

Proof.
Assume i e   with be a root of 0 π  .Substituting this root to Equation (11) yields Separating the real and imaginary parts, we have and therefore we get  Proof.From Equations ( 11) and ( 15) we have Based on Theorem 1 and Lemmas 1-3, we have the following results on stability and bifurcation of system (5) at the DFE.

1) If and
then the DFE is asymptotically stable for all 0   .Proof.

1) Let and
then it is known from Lemmas 1 and 2 that the characteristic Equation ( 10) has no root with modulus one for all 0   .Applying Theorem 1, all roots of Equation ( 10) have modulus less than one for all 0   .Thus, conclusion (a) follows.
2) Let and then one of roots of characteristics Equation (10) has modulus less than one for all Other roots of Equation (10) satisfy Equation (11).From Lemmas 1 and 3 we know that all roots of Equation (11) have modulus less than one when  0 0,     , and Equation (11) has at least a couple of roots with modulus greater than one when 0    .Hence conclusion (b) fol-lows.

Endemic Equilibrium
The linearization of system (7)  , where 1 ln Using the same transformation as in the DFE, i.e. 16) can be written as where . The constant matrix B in Equation ( 17) is It is easy to show that the characteristic equation of matrix B is where   Clearly this characteristic equation has a trivial root and other roots are determined by .Direct calculations show that if is less than one for all 0   and other roots satisfy where   and   satisfy Equation (1 .nd 5) Theorem 3. s 2.1 and 3.1 give the sufficient conditions for DFE and EE to be asymptotically stable, respectively.iscrete sy em obtained by the Eule Theorem Unlike the d st r method where the stability conditions for equilibria depend on the size of time-step h, the proposed discrete system (7) has stability conditions which are consistent with the stability conditions for equilibria of continuous system (4) for any size of time step h; see [7] for the stability properties of the continuous system.In addition, Theorems 2.2 and 3.2 show the existence of bifurcation controlled by the time delay  .

Numerical Simulations
To confirm our pre io v us theoretical analysis, in this secsimulations using non-) with time -step h = tion we present some numerical standard finite difference scheme (7 0.1 (or k = 10).For the simulations we adopt parameters used by [7], i.e.

Endemic Equilibrium (R 0 > 1)
For the EE sce e such that R 0 = the numerical nario we us 5.0   2.3810 > 1. Figures 3(a then the behavior of numerical solutions changes from a solution converging to the EE (0.7577, 1.3064) to a solution which oscillates periodically about the EE.Finally we compare the numerical solutions obtained by Euler method with those obtained by our NSFD method.We found that Euler method with a relatively big numerical time-step (h) may produce unrealistic negative and unbounded solutions.However, using the same or even bigger numerical time-step, NSFD method always gives positive and bounded solution.Depending on the parameters used in the simulations, the solution will be periodic or convergent to a correct equilibrium point.For example, in Figure 5 we plot the results of Euler method and NSFD method for the EE scenario using h = 0.25.It is seen from

Conclusion
In this paper we have introduc epidemic model o ference method.U scheme reproduces exactly the same equilibria as well as their stability conditions as those of continuous model, i.e. if are the death rates of immature and mature population, respectively.The delay  is considered as the maturation time.The parameter 0 respectively the constant contact rate, the disease induced death rate constant and the recovery constant rate, respectively.In this paper we consider a special class of Equation (1) by assuming that the death rate in each stage prior to the adult stage and the disease induced death are negligible, i.e. 1 0     .Using the Rick function as the birth rate function


. This root depends continuously on  .Based on Equation(11) we have where   and   satisfy Equation(15).

Figure 3 .
Figure 3.The numerical solutions of the discrete SIS model with a delay for the EE scenario with 1 < p/δ < e 2 and R 0 > 1 using (a) τ = 1.0 and (b) τ = 4.0.

Figure 4 .
Figure 4.The numerical solutions of the discrete SIS model with a delay for the EE scenario with p/δ > e 2 , R 0 > 1 using (a) τ = 1.60 and (b) τ = 1.70.

Figure 5 .
Figure 5.The numerical solutions for the EE scenario with p/δ > e 2 , R 0 > 1, h = 0.25 (k = 4) and τ = 2 obtained by (a) Euler method and (b) NSFD method.Note that the result o does not affect the stability of the equilibria.However, in the case of 2 p e   , the stability of equilibria changes when the delay p al value.Here the discrete SIS model with a delay has periodic solutions when the stability is her words, a Neimark-Sacker bifurcation occurs.asses a critic lost.In ot around the EE  