Dynamics of a SIR Epidemic Model of Childhood Diseases with a Saturated Incidence Rate: Continuous Model and Its Nonstandard Finite Difference Discretization

A SIR epidemic model that describes the dynamics of childhood disease with a saturated incidence rate and vaccination program at a constant rate was investigated. For the continuous model we first show its basic properties, namely, the non-negativity and boundedness of solutions. Then we investigate the existence and both local and global stability of the equilibrium points. It was found that the existence and stability properties of equilibrium points fully determined the basic reproduction number. We also propose and analyze a discrete-time analogue of the continuous childhood diseases by applying a nonstandard finite difference method. It is shown that our discrete model preserves the dynamical properties of the corresponding continuous model, such as the positivity solutions, the population conservation law, the existence of equilibrium points and their global stability properties.


Introduction
One of the major issues in public health problems is childhood diseases (measles, mumps, influenza, smallpox, chickenpox, Rubella, Polio, etc.), which can spread rapidly among children age 5 and below due to frequent contact with their peers at school or elsewhere. For childhood diseases, vaccination is considered as an effective control strategy against childhood diseases. To understand the transmission dynamics of childhood diseases, Makinde [1] considered a SIR epidemic model which includes three sub-populations, namely, the susceptible sub-population (S), the infective sub-population (I) and the recovered sub-population (R). The influence of a preventive vaccine is modeled by assuming that there exists a fraction of the susceptible population vaccinated at birth each year. The incidence rate of infection in this model is assumed to be bi-linear, i.e., bSI, where b is a parameter that measures the disease transmission rate. However, when an epidemic occurs in a large population, there are changes in people's behavior to prevent infection (known as the crowding effect [2][3][4][5]) and hence the bi-linear incidence rate cannot be applied anymore. To consider such behavioral changes, Anderson and May [6] proposed a saturated incidence rate 1 1+aI , where a is the inhibitory constant. The infection rate clearly decreases as the constant of inhibition increases, and in response to increasing the preventive measures and proper awareness.
In this paper, we modify the model in [1] by considering a saturated incidence rate. The compartment diagram of the proposed model is illustrated in Figure 1. One of the important preventive measures in the spread of childhood disease is vaccination. Therefore, we introduce an additional parameter θ to represent the fraction of individuals vaccinated at birth each year where 0 ≤ θ ≤ 1. The SIR epidemic model of childhood diseases with saturated incidence rate and constant vaccination strategy is where A is the constant recruitment rate; m and e are the natural death rate and the death rate caused by the disease, respectively; and g represents the recovery rate. It is assumed that the initial values of system (1) are non-negative, namely, S(0) = S 0 ≥ 0, I(0) = I 0 ≥ 0 and R(0) = R 0 ≥ 0. From a mathematical point of view, Makinde [1] has provided a qualitative analysis for the SIR model with a bi-linear incidence rate. He also applied the Adomian decomposition method to obtain an analytical approximation to its solutions. Another analytical approximation to this problem has also been introduced by Yildirim and Cherruault [7], namely, by using the homotopy perturbation method. On the other hand, it is often desirable to get quantitative solutions and therefore the epidemic model has to be discretized. In this case, the obtained numerical scheme has to maintain the dynamical properties of the corresponding continuous model. For this purpose, Mickens [8] has introduced a nonstandard finite difference (NSFD) method, which has two characteristics: The first order derivative in the continuous model is approximated by the generalized forward difference The nonlinear terms are approximated non-locally.
For all we know, the dynamics of the SIR epidemic model (1) and its discrete version have not been studied. Hence, the SIR epidemic model of childhood diseases with saturated incidence rate and constant vaccination strategy is proposed and the dynamics are investigated in this paper. We first discuss the basic properties of the model such as the non-negativity solutions and the population conservation law. We also provide local and global stability analysis for all possible equilibrium points. Then, we construct an NSFD scheme for the model (1) and investigate the dynamics of the obtained discrete model. The dynamics of the discrete model are compared to the continuous model to check their consistency.

Basic Properties of the Continuous Model
Since we are considering a model of population dynamics, all solutions of model (1) must be non-negative. Furthermore, due to the limited resources, the size of the population has to be bounded. In the following theorem, we prove the non-negativity and the boundedness of solutions of system (1). Theorem 1. All solutions of system (1) with non-negative initial values are always non-negative and uniformly bounded. Proof. From system (1) and the non-negativity of initial values S 0 ≥ 0, I 0 ≥ 0 and R 0 ≥ 0, we have Hence all solutions of system (1) are always non-negative. Next, if we define the total population as T(t) = S(t) + I(t) + R(t), then system (1) gives the following population conservation law (see [14,21] for the detail terminology of population conservation law): Using the Gronwall's inequality, we get where T 0 = S 0 + I 0 + R 0 . Hence lim t→+∞ T(t) ≤ A m , which shows that the solutions of system (1) are uniformly bounded.
It is easy to show that system (1) always has a disease free equilibrium (DFE) point E 0 = S 0 , 0, R 0 , where S 0 = (1 − θ)A/m and R 0 = θ A/m. Furthermore, if R 0 > 1, then the system also has a unique endemic equilibrium (EE) point E * = (S * , denotes the basic reproduction number which can be determined easily using the next generation matrix method; see [22]. In the following part, we investigate the local and global stability of the existing equilibrium points. Since the first two equations of system (1) do not depend on R, it is enough to consider the following reduced system

Local Stability Analysis of the Continuous Model
To study the local stability of each equilibrium point, we examine the eigenvalues of the Jacobian matrix of the system (6) at each equilibrium point. The Jacobian matrix of system (6) at E 0 = (S 0 , 0) is The eigenvalues of J(E 0 ) are λ 1 = −m < 0 and λ 2 = bS 0 − (m + g + e). We observe that λ 2 < 0 if R 0 < 1 and λ 2 > 0 whenever R 0 > 1. Hence, the DFE is locally asymptotically stable (L.A.S.) if R 0 < 1 and it is unstable if R 0 > 1. Next, we note that the Jacobian matrix at the EE point E * = (S * , I * ) is given by from which we get Thus, the real parts of all eigenvalues of the system(6) are negative, and consequently, the EE point is L.A.S. The local stability of equilibrium points of system (6) is summarized in the following theorem.

Global Stability Analysis of the Continuous Model
In this section, we discuss the global stability of equilibrium points. The global stability of the DFE point is shown by choosing a suitable Lyapunov function and employing the La Salle invariance principle, while that of the EE point is proven by defining a Dulac function and applying the Poincaré-Bendixson theory. The sufficient conditions for the global stability of both equilibrium points are stated in the following theorems. (1), the time derivative of our Lyapunov function is given by Proof. It has been shown previously that E * is L.A.S. if R 0 > 1. Next, we show the system (6) has no closed-orbit. For that, we define a Dulac function as follows: Then, we have that Thus, system (6) or equivalently system (1) does not have a closed orbit. Since the EE point

Construction of the Discrete NSFD Model
To get a dynamically-consistent numerical scheme, in this section we apply the NSFD method to discretize system (1). Here, the numerical approximations of S(t), I(t) and R(t) at t = nh, n = 0, 1, 2, . . . are respectively denoted by S n , I n and R n , where h is the time-step size. Following Mickens [8], we employ a generalized forward difference and nonlocal approximation for system (1): where ϕ(h) is the denominator function. As in the previous continuous model, we also consider a non-negative initial values S 0 ≥ 0, I 0 ≥ 0 and R 0 ≥ 0. By adding all equations in the scheme (9), we get where T n = S n + I n + R n . It can be seen that Equation (10) is the approximation for the continuous population conservation law (3). To satisfy the exact population conservation law, we choose the following denominator function: Notice that we have implemented a nonlocal approximation for the right hand sides of system (1) and the denominator function satisfies ϕ(h) = h + O(h 2 ) > 0 for all h > 0. Hence, scheme (10) and scheme (9) are considered to be the NSFD scheme. Furthermore, if the denominator function (11) is substituted into the inequality (10) then we can easily show that Based on inequality (12), we can prove using mathematical induction that for any h > 0 the solution of Equation (10) with the denominator function (11) satisfies the exact population conservation law (4): The NSFD scheme (9) is implicit, but it can be written in an explicit form as follows: where Φ(x) = bx 1+ax and Φ n = Φ(I n ). By remembering that the initial values are non-negative, ϕ(h) > 0 and all parameters are positive with 0 ≤ θ ≤ 1, it is easily seen that Φ 0 ≥ 0 and S 1 ≥ 0. Consequently, we also have I 1 ≥ 0 and R 1 ≥ 0. By applying this procedure consecutively for all n, we conclude that numerical solutions of system (14) are always non-negative: S n ≥ 0, I n ≥ 0 and R n ≥ 0 for all n > 0 and for any h > 0.
Using standard calculations, we can verify that the discrete NSFD model (14) always has a DFE point E 0 = (S 0 , 0, R 0 ). If R 0 > 1, then the discrete system (14) also has a unique EE point E * = (S * , I * , R * ). These equilibrium points and their existence conditions are exactly the same as those of the continuous model (1), irrespective of h. In the following sub-sections we investigate the stability of those equilibrium points.

Local Stability Analysis of the Discrete NSFD Model
We notice that the first two equations in system (14) are independent of variabel R n . Hence, for the stability analysis, we can consider the following reduced discrete NSFD model .
The DFE and EE points of the reduced discrete NSFD model (15) can respectively be written as E 0 = (S 0 , 0) and E * = (S * , I * ). For simplification of analysis, we introduce the following functions: The Jacobian matrix of system (15) at an equilibrium pointÊ = (Ŝ,Î) is The local stability of the DFE (E 0 ) and the EE point are stated in the following theorems.

Proof. The Jacobian matrix (16) at E 0 is
It is clear that the eigenvalues of J(E 0 ) are Using the definition of R 0 , see (5), we can show that if R 0 < 1 then 0 < λ 2 < 1, and thus E 0 is L.A.S. On the contrary, it is obvious to verify that λ 2 > 1 if R 0 > 1, which shows that E 0 is unstable. Notice that the stability properties of E 0 are independent of h. Proof. The Jacobian matrix (16) where . Clearly, p > 1, q > 1, r > 0, s > 0 and p > r. The characteristic equation of J(E * ) is given by where P = 1 p + 1 p 2 q (p 2 + s(p − r)) > 0 and Q = p+s p 2 q > 0. By performing some straightforward calculations, we get the following results.

•
Since p > 1, we have From the definition of p, q, r and s, we have , we directly get Based on the Schur-Cohn criterion [22,23], we conclude that the roots of characteristic Equation (19) are inside the unit disk and therefore E * is L.A.S. if R 0 > 1 for any h.

Global Stability Analysis of the Discrete NSFD Model
We now investigate the global stability of each equilibrium point using a suitable Lyapunov function. Theorem 7. If R 0 ≤ 1, then the DFE point (E 0 ) of the discrete NSFD model (15) is GAS.
Proof. Let us consider a discrete Lyapunov function where f (z) = z − 1 − ln(z) ≥ f (1) = 0. Using the fact that ln(z) ≤ z − 1, we get As R 0 ≤ 1, it is clear that U n+1 − U n ≤ 0 for all n ≥ 0. This shows that U n is monotonic decreasing sequence. By noting that U n ≥ 0, we have that lim n→∞ U n ≥ 0 exists and lim n→∞ (U n+1 − U n ) = 0. Thus, we get lim n→∞ S n+1 = S 0 and lim n→∞ (1 − R 0 )I n+1 = 0. It is obvious that if R 0 < 1 then lim n→∞ I n+1 = 0. Moreover, if R 0 = 1, then according to the first equation of the discrete NSFD model (15) and using the fact that lim n→∞ S n+1 = S 0 , we get lim n→∞ I n = 0. Then we conclude that E 0 is GAS.
It has been mentioned that if R 0 > 1 then, besides the DFE point, the reduced discrete NSFD model (15) has also an EE point E * . Next, we prove that whenever E * exists then it is GAS.
Proof. To prove the global stability of point (E * ), we follow the trick of Enatsu et al. [24], namely, by considering a discrete Lyapunov function where V 1 n = S * f S n S * , V 2 n = I * f I n I * , V 3 n = f Φ n Φ * and f (z) is defined as in the proof of Theorem 7. Here we denote Φ * = Φ(I * ). Based on the definition of V 1 n and V 2 n , we have and From (22) and (23), we get Obviously that V n is monotonic decreasing sequence since V n+1 − V n ≤ 0 for all n ≥ 0. As V n ≥ 0, lim n→∞ V n ≥ 0 exists and consequently lim n→∞ (V n+1 − V n ) = 0. Then we have lim n→∞ S n+1 = S * and lim n→∞ I n+1 = I * . Hence, the EE point E * = (S * , I * ) is GAS.

Numerical Simulations
To illustrate our previous analytical results, we provide some results of numerical simulations using hypothetical parameters:  Finally, we performed numerical simulations using the same parameter values as in Figure 3 but with different initial values, namely, S(0) = 10, I(0) = 6 and R(0) = 1. We plotted numerical solutions obtained by the proposed NSFD scheme (9) and the explicit Euler scheme in Figure 4. This figure shows that the NSFD scheme gives non-negative solutions that are convergent to the correct equilibrium point. Although the solutions obtained by the explicit Euler scheme with h = 0.5 and h = 0.75 are also convergent to the correct equilibrium point, the solutions are biologically unrealistic because the solutions are negative at some points. The explicit Euler scheme with higher time-step sizes may also give a non-convergent solution; see Figure 4h. Figures 3 and 4 confirm that the proposed NSFD scheme (9) is dynamically consistent with the model (1) regardless the value of h. On the contrary, the dynamical consistency of the explicit Euler scheme was achieved only for a relatively small h.

Conclusions
A SIR epidemic model of childhood diseases with a saturated incidence rate and constant vaccination strategy has been introduced. The non-negativity of the solutions and the conservation of population law for the proposed model have been shown. It is also shown that the model has always a DFE point which is (locally and globally) asymptotically stable if R 0 < 1. If R 0 > 1, then the DFE point is unstable and there appears an EE point, which is (locally and globally) asymptotically stable. In this paper, we have also proposed a discrete-time analog of the continuous childhood diseases epidemic model. To achieve this objective, we discretized the continuous model using the NSFD method. It is shown that the solutions of the constructed discrete NSFD model are always non-negative and satisfy the exact population conservation law. Furthermore, by verifying the eigenvalues of the Jacobian matrix of the proposed discrete NSFD model and taking suitable discrete Lyapunov functions, we have established the local and global stability for both DFE and EE points. Both the local and stability properties of the discrete NSFD model are exactly the same as those of the corresponding continuous model, which are entirely controlled by the basic reproduction number and do not depend on the time-step size. Hence, the proposed NSFD discrete model is dynamically consistent with the continuous SIR epidemic model of childhood diseases. The analytical results have been confirmed by our numerical simulations.