On the nonstandard numerical discretization of SIR epidemic model with a saturated incidence rate and vaccination

1.   Introduction

Over recent decades, mathematical models are considered as crucial tools for understanding the transmission mechanism and evaluating the control strategies of infectious diseases, see ^[1,2,3]. In this respect, Yusuf and Benyah ^[4] have introduced and analyzed a SIR epidemic model with constant vaccination and treatment where the incidence rate is assumed to be bi-linear. The bi-linear incidence rate describes that disease transmission is proportional to the number of susceptible individuals and infective individuals, and so the greater the number of susceptible individuals and infective individuals cause an even greater transmission rate. Hence the bi-linear incidence rate is not realistic for cases in a large population since there are changes in the social and psychological behavior of society or some prevention measure taken by the society which reduce the transmission rates. To deal with the crowding effect or behavioral change, Capasso and Serio ^[5] proposed a saturated incidence rate $\frac{\beta S I}{1+\alpha I}$ where $\beta$ is the maximal disease transmission rate and $\alpha$ is the saturation factor that measures the inhibitory effect due to the crowding effect or behavioral change of society ^[6]. After that, the saturated incidence rate has been widely applied to deterministic epidemic models ^[7,8,9,10], epidemic models with optimal control ^[11,12,13] or stochastic epidemic models ^[14,15,16]. If the bilinear incidence rate in the SIR epidemic model ^[4] is replaced by the saturated incidence rate, then we get

where $S \equiv S(t), I \equiv I(t)$ and $R \equiv R(t)$ denote the sub-population of susceptible, infective, and recovered at time $t$, respectively. $\Lambda, \nu$ and $\gamma$ are positive parameters which respectively represent the recruitment rate of population, vaccination rate, and recovery rate due to treatment. The natural death rate and the death rate induced by the disease are respectively denoted by $\mu$ and $\delta$. It is assumed that initial values are non-negativity: $S(0) = S_0 \ge 0, I(0) = I_0 \ge 0$ and $R(0) = R_0 \ge 0$. Recently, the SIR model (1.1) has been applied by Khan et al. ^[17] to model the spread of hepatitis B virus. The local and global stability analysis of the disease free equilibrium (DFE) point were established, but the disease endemic equilibrium (DEE) point was only studied locally. The global stability analysis of the DEE point of the model was then provided by Hoang and Egbelowo ^[18]. Very recently, Macìas-Dìaz et al. ^[19] studied the diffusion effect on the spread of the hepatitis B virus.

The mathematical epidemic models are generally in the form of a nonlinear differential equations system so that their analytical solutions are difficult to be obtained. Hence numerical approximations are often needed to get reliable results. However, it was shown that many standard numerical methods such as the Euler method, Runge-Kutta method, and other standard finite difference methods may fail to maintain the dynamical properties of the corresponding continuous models, see e.g. ^[20,21,22]. To overcome such dynamics-inconsistency, Mickens ^[23] has developed a nonstandard finite difference (NSFD) method. A numerical scheme for a system of first-order differential equations is called an NSFD scheme if at least one of the following conditions is satisfied.

● The first-order derivatives in the system are approximated by the generalized forward difference method

where $u_n = u(t_n)$ and $\phi \equiv \phi(h)$ is the denominator function such that $\phi(h) = h + \mathcal{O}(h^2)$ and $h$ is the time-step size.

● The nonlinear terms are approximated non-locally. For example $u^2(t_n) \approx u_n u_{n+1}$.

The NSFD method has been successfully implemented in various problems such as epidemiology ^{[24,25,26,27,28,29,30]}, ecology ^{[28,31,32,33]} or metapopulation model ^[34]. However, we notice that an NSFD scheme for a differential equation is not unique. Hence, the dynamical properties of an NSFD scheme need to be analyzed to ensure whether an NSFD scheme maintains the dynamics-consistency with its original continuous equation or not. NSFD schemes in ^{[27,28,31,32,33]} did not apply the generalized forward difference method but they used nonlocal approximation. Moreover, the authors in those references did not consider the population conservation law. Their dynamical analysis shows that the NSFD schemes in ^[27,31,32] maintain the stability properties of all existing equilibrium points irrespective of the time-step size ($h$). However, the NSFD schemes in ^[28,33] require a relatively small $h$ for the stability of equilibrium points. Cui et al. ^[24,25,26] and the authors in ^[29,30] applied both generalized forward difference method and nonlocal approximation. They have proven that the obtained NSFD schemes preserve all dynamical properties of the corresponding continuous models irrespective of $h$. In particular, the denominator function must be appropriately chosen to satisfy the exact population conservation law.

Recently, Hoang and Egbelowo ^[18] proposed a numerical scheme for the model (1.1) to get the following discrete model

where $\phi$ is any denominator function but it must obey $\phi \equiv \phi(h) = h + \mathcal{O}(h^2)$. Notice that the NSFD scheme (1.2) implements the generalized forward method and a local approximation for the right-hand sides. Rigorous investigation on the dynamical properties shows that the discrete model (1.2) preserves the dynamics-consistency of the model (1.1) only if the denominator function satisfies $\phi(h) < \phi^{\ast}$, where $\phi^{\ast}$ is a critical value of the denominator function which is determined by some parameters in the model (1.1). Hence, the preservation of dynamics-consistency is achieved only when the time-step size is relatively small.

In this paper, we propose an alternative NSFD scheme for the model (1.1) which is dynamically-consistent with the original continuous model for any time-step size. Following Cui et al. ^{[24,25,26,29,30]}, we apply both the generalized forward difference method and nonlocal approximation to construct the NSFD scheme. We organize the rest of the paper as follows. In Section 2 we review the dynamical properties of model (1.1) based on the analytical results in ^[17,18]. The derivation of the NSFD scheme for model (1.1) is presented in Section 3. We show in this section that the obtained discrete NSFD model has the same dynamical properties as the corresponding continuous model. To illustrate our analytical findings, we present numerical simulations in Section 4. We draw some conclusions in the final section.

2.   Dynamics of the continuous epidemic model

As for other models of population dynamics, the SIR epidemic model with a saturated incidence rate (1.1) is required to have positive and bounded solutions. The positivity of the solution of model (1.1) has been proven in ^[17] and is stated in the following theorem.

Theorem 1. ^[17] If the solution of the SIR epidemic model (1.1) subject to non-negative initial values exist, then it is positive $\forall t > 0$.

The boundedness of the solution is related to limited resources so that the population cannot grow infinitely. Hence, the boundedness of the solution of model (1.1) is biologically important but it has not been described in ^[17,18]. In the following theorem, we establish the boundedness of solution for the model (1.1).

Theorem 2. The solution of the SIR epidemic model (1.1) subject to non-negative initial values is ultimately bounded.

Proof. Let $P(t) = S(t)+I(t)+R(t)$ be the total population. From model (1.1) we have

from which we get

where $P_0 = S_0+I_0+R_0$. It is clear that $\lim_{t \to +\infty} P(t) \le \frac{\Lambda}{\mu}$, and thus $S(t), I(t)$ and $R(t)$ are ultimately bounded.

From Theorem 1 and Theorem 2, we see that model (1.1) has a feasible region

which is positively invariant. Therefore, the SIR epidemic model (1.1) is mathematically and epidemiologically well posed in $\Omega$.

Using the next generation matrix method, Khan et al. ^[17] derived the basic reproduction number for model (1.1), namely

It was also shown in ^[17] that model (1.1) always has a DFE point $E^0 = \left(S^0, 0, R^0\right)$, where $S^0 = \frac{\Lambda}{\mu+\nu}$ and $R^0 = \frac{\Lambda\nu}{\mu(\mu+\nu)}$. In addition to the DFE point, if $\mathcal{R}_0 > 1$, then model (1.1) also has a DEE point $E^{\ast} = \left(S^{\ast}, I^{\ast}, R^{\ast}\right)$ where $S^{\ast} = \frac{\alpha \Lambda + (\mu+\gamma+\delta)}{\alpha (\mu+\nu)+\beta}, I^{\ast} = \frac{(\mu+\nu) (\mathcal{R}_0 - 1)}{\alpha (\mu+\nu)+\beta}$ and $R^{\ast} = \frac{\nu S^{\ast} + \gamma I^{\ast}}{\mu}$. The following theorem about the stability of the DFE and DEE points are taken from ^[17,18].

Theorem 3. ^[17]

1. The DFE point of model (1.1) is locally asymptotically stable if $\mathcal{R}_0 < 1$ and it is unstable if $\mathcal{R}_0 > 1$.

2. The DEE point of model (1.1) is locally asymptotically stable if $\mathcal{R}_0 > 1$.

Theorem 4. ^[18]

1. The DFE point of model (1.1) is globally asymptotically stable if $\mathcal{R}_0 \le 1$.

2. The DEE point of model (1.1) is globally asymptotically stable if $\mathcal{R}_0 > 1$.

3.   Dynamically-consistent NSFD discretization of the SIR epidemic model with saturated incidence rate and vaccination

3.1.   Construction of the NSFD scheme

For the rest of this paper, we respectively denote $S_n, I_n$ and $R_n$ as the numerical approximation of $S(t), I(t)$ and $R(t)$ at $t = nh, n = 0, 1, 2, ...$, where $h$ is the time-step size. Using the idea of Mickens ^[23], we discretize model (1.1) as follows

The initial values of the discrete NSFD SIR epidemic model (3.1) are also supposed to be non-negative: $S_0 \ge 0, I_0 \ge 0$ and $R_0 \ge 0$. Different from the scheme (1.2), the NSFD scheme (3.1) approximates the nonlinear terms in the right-hand sides nonlocally. Furthermore, the denominator function $\phi$ in (3.1) is not arbitrary, but it will be derived such that it obeys the exact boundedness of solutions. To do so, we denote the total population as $P_n = S_n+I_n+R_n$. From the discrete NSFD model (3.1), we have

By straightforward calculations, it is easy to verify that if we take a denominator function

then the solution of equation (3.2) satisfies

for any $h > 0$. Additionally, for the case of $\delta = 0$, we can also show that the total population of the discrete NSFD model (3.1) are exactly the same as that of model (1.1). Hence, the discrete NSFD model (3.1) maintains the exact boundedness of solutions (2.2). The discrete NSFD model (3.1) can be rearranged to get explicit form

where $\Psi(z) = \frac{\beta z}{1+\alpha z}$ and $\Psi_n = \Psi(I_n)$. Since all parameters in (3.5) are positive, it is proven that $S_n \ge 0, I_n \ge 0$ and $R_n \ge 0$ for all $n > 0$ and for any $h > 0$. This finding confirms the discrete NSFD model (3.5) preserves the non-negativity solutions for any $h$.

By applying some algebraic calculations, we can easily verify that the discrete NSFD model (3.5) has exactly the same equilibrium points as continuous model (1.1), namely the DFE point $E^0 = (S^0, 0, R^0)$ and the DEE point $E^{\ast} = \left(S^{\ast}, I^{\ast}, R^{\ast}\right)$. The DEE point $E^{\ast}$ exists if $\mathcal{R}_0 > 1$. The local and global stability of equilibrium point $E^0$ and $E^{\ast}$ are investigated in the following Subsections.

3.2.   Local stability analysis of the discrete NSFD model

By noticing that the first two equations in the discrete NSFD model (3.5) do not depend on the third equation, the stability analysis can be performed by considering the following reduced model

The Jacobian matrix of the discrete NSFD model (3.6) at a point $\bar{E} = (\bar{S}, \bar{I})$ is given by

By observing the eigenvalues of the Jacobian matrix $J(\bar{S}, \bar{I})$, we study the local stability of equilibrium points as in the following theorems.

Theorem 5. If $\mathcal{R}_0 < 1$, then the DFE point $E^0$ of the discrete NSFD model (3.6) is locally asymptotically stable for all $h > 0$. Furthermore, if $\mathcal{R}_0 > 1$ then $E^0$ is unstable for all $h > 0$.

Proof. If we substitute the DFE point $E^0$ into the Jacobian matrix (3.7), the we get

The eigenvalues of $J(E^0)$ are obviously

Clearly that if $\mathcal{R}_0 < 1$ then $0 < \lambda_2 < 1$, while if $\mathcal{R}_0 > 1$ then $\lambda_2 > 1$. This fact confirms that if $\mathcal{R}_0 < 1$ then the DFE $E^0$ is locally asymptotically stable and it is unstable if $\mathcal{R}_0 > 1$, irrespective of $h$.

To prove the local stability of DEE $E^{\ast}$, we need the following Schur-Cohn criterion ^[1,35].

Lemma 6. ^[1,35] The roots of quadratic equation $\lambda^2 - T \lambda + D = 0$ satisfy $\left|\lambda_i\right| < 1, i = 1, 2$ if and only if all of the following conditions are fulfilled:

1. $D < 1$,

2. $1+T+D > 0$,

3. $1-T+D > 0$.

Theorem 7. If $\mathcal{R}_0 > 1$, then the DEE point $E^{\ast}$ of the discrete NSFD model (3.6) is locally asymptotically stable for all $h > 0$.

Proof. The Jacobian matrix (3.7) evaluated at the DEE point $E^{\ast}$ is

where $a_1 = 1+\phi(\Psi^{\ast}+\mu+\nu) > 1, \; \; a_2 = 1+\phi(\mu+\gamma+\delta) > 1, \; \; a_3 = \phi\Psi^{\ast} > 0$ and $\; \; a_4 = \frac{\phi \beta (S^{\ast}+\phi\Lambda)}{(1+ \alpha I^{\ast})^2} > 0$. The eigenvalues of $J(E^{\ast})$ is determined by the following quadratic equation

where $T = \text{Trace}\left(J(E^{\ast})\right) = \frac{1}{a_1} + \frac{1}{a_1^2 a_2}(a_1^2+a_4(a_1-a_3))$ and $D = \text{Determinant}\left(J(E^{\ast})\right) = \frac{a_1+a_4}{a_1^2 a_2}$. Since $a_1 > a_3$, it is apparently that $T > 0$ and $D > 0$. We also have the following results.

1. Due to $a_1 > 1$, we get

2. It is obvious that $1+T+D > 0$.

3. Straightforward calculations show that if $\mathcal{R}_0 > 1$ then

Hence, the conditions for the Schur-Cohn criterion in Lemma 6 are satisfied whenever $\mathcal{R}_0 > 1$. It is concluded that the DEE point $E^{\ast}$ is locally asymptotically stable if $\mathcal{R}_0 > 1$ for any $h$.

3.3.   Global stability analysis of the discrete NSFD model

This Subsection is devoted to establish the sufficient conditions for the global stability of equilibrium points of the discrete NSFD model (3.6).

Theorem 8. If $\mathcal{R}_0 \le 1$, then the disease free equilibrium point $E^0$ of the discrete NSFD model (3.6) is globally asymtotically stable for all $h > 0$.

Proof. We first introduce a discrete Lyapunov function

where $\varphi(z) = z - 1 - \ln(z) \ge \varphi(1) = 0$. Based on the first two equations of (3.1) we obtain

Obviously that if $\mathcal{R}_0 \le 1$, then $\mathcal{U}_{n+1} - \mathcal{U}_n \le 0$ for all $n \ge 0$, which shows $\mathcal{U}_n$ is monotone decreasing sequence. Since $\mathcal{U}_n \ge 0$, there exists a $\lim_{n\to\infty} \mathcal{U}_n \ge 0$ and thus $\lim_{n\to\infty} \left(\mathcal{U}_{n+1} - \mathcal{U}_n\right) = 0$. Consequently we have $\lim_{n\to\infty} S_{n+1} = S^0$ and $\lim_{n\to\infty} \left(1 - \frac{\mathcal{R}_0}{1+\alpha I_{n+1}} \right) I_{n+1} = 0$ for any $h$. It is straightforward to show that if $\mathcal{R}_0 \le 1$ then $\lim_{n\to\infty} I_{n+1} = 0$. This completes the proof.

In the following theorem we show that $E^{\ast}$ is globally asymtotically stable if it exists.

Theorem 9. If $\mathcal{R}_0 > 1$, then the DEE point $E^{\ast}$ of the discrete NSFD model (3.6) is globally asymptotically stable for all $h > 0$.

Proof. To establish the sufficient condition for which $E^{\ast}$ is globally stable, we apply the approach of Enatsu et al. ^[36]. Here we define the following discrete Lyapunov function

where $\Psi^{\ast} = \Psi(I^{\ast})$ and $\varphi(z)$ is defined as in the proof of Theorem 8. By substituting the first two equations in (3.1) into equation (3.12), we obtain

By denoting $\chi_n = \frac{S_n}{S^{\ast}}, \xi_n = \frac{I_n}{I^{\ast}}$ and $\zeta_n = \frac{\Psi_n}{\Psi^{\ast}}$, equation (3.13) can be written as

By definition of $\varphi(z)$, we can verify that

Therefore $\mathcal{W}_n$ is monotone decreasing sequence. Using the same argument as in the proof of Theorem 8 we can show that $\lim_{n\to\infty} \left(\mathcal{W}_{n+1} - \mathcal{W}_n\right) = 0$, from which we get $\lim_{n\to\infty} S_{n+1} = S^{\ast}$ and $\lim_{n\to\infty} I_{n+1} = I^{\ast}$ for all $h$. Hence, we have the theorem.

4.   Numerical simulations

To give better view of the previous analytical results, we present some numerical simulations using the proposed discrete NSFD model (3.5). For the simulations, we take initial value $S(0) = 4, I(0) = 0.2$, and hypothetical parameters: $\Lambda = 4, \mu = 0.5, \beta = 0.5, \alpha = 0.01, \gamma = 0.1, \nu = 0.2$, and $\delta = 0.1$. By calculation, the basic reproduction number is $\mathcal{R}_0 = 4.0816 > 1$. The equilibrium points of both continuous model (1.1) and discrete NSFD model (3.5) are exactly the same, namely the unstable DFE point $E^0 = (5.7143, 0.0, 2.2857)$ and the globally asymptotically stable DEE point $E^{\ast} = (1.4590, 4.2547, 1.4348)$. In Figure 1 we compare the solutions of the proposed discrete NSFD model (3.5) with those of the explicit Euler method using $h = 1$ and $h = 1.35$. It is shown that the explicit Euler method produces a periodic solution of period two for the case of $h = 1$ and a chaotic behavior for $h = 1.35$. In contrast to the explicit Euler method, the solutions of the proposed NSFD scheme (3.5) are convergent to the correct equilibrium point $E^{\ast}$ for both $h = 1$ and $h = 1.35$. The complete dynamical behavior of the numerical results of the proposed discrete NSFD model (3.5) and the explicit Euler method with respect to $h$ can be seen in Figure 2. Figure 2.(a-b) indicates that there exists a critical time-step size $h^{\ast}$ such that the solution obtained by the explicit Euler method is convergent to the stable DEE point $E^{\ast}$ if $h < h^{\ast}$. For larger time-step size, the solution may converge to a stable periodic solution of period two or that of period four, or even exhibits a chaotic behavior. On the contrary, the solutions obtained by the proposed NSFD scheme (3.5) are always convergent to the DEE point $E^{\ast}$ for any $h$, see Figure 2.(c-d). These numerical simulations confirm that the dynamics of the discrete model obtained by the explicit Euler method are dependent on the time-step size, while the dynamics of the proposed discrete NSFD model (3.5) are independent of $h$.

Next, using the same parameter values as before, we compare the numerical solutions obtained by the NSFD scheme (1.2) ^[18] and our proposed NSFD scheme (3.5), see Figure 3. Here, the denominator function for the NSFD scheme (1.2) is $\phi = \frac{1-\exp(\mu h)}{\mu}$. As stated in the Introduction, the NSFD scheme proposed in ^[18] is dynamically consistent only for restrictive time-step size. Indeed, according to the analytical results in ^[18], the NSFD scheme (1.2) is dynamically-consistent if $\phi^{\ast} < 0.9486$ which is equivalent to $h < 1.2860$. Figure 3 confirms this phenomenon. If we take $h = 1.5$ then the numerical solution obtained by the NSFD scheme (1.2) does not converge to the DEE point but it goes to a periodic solution of period two, see Figure 3.(a). If we further increase the time-step size $(h = 2)$ then we get a chaotic solution, see Figure 3.(b). On the other hand, the proposed NSFD scheme (3.5) preserves the stability of the DEE point, even for relatively large $h$, see Figure 3.(c-d).

Finally, we perform simulation using parameters $\Lambda = 1, \mu = 0.5, \beta = 0.25, \alpha = 0.01, \gamma = 0.1, \nu = 0.2, \delta = 0.1$, and time step size $h = 0.1$ with initial value $S(0) = 6$ and $I(0) = 1$. The reproduction number in this case is $\mathcal{R}_0 = 0.5102 < 1$. Hence, the DFE point $E^0 = (1.4286, 0.0, 0.5714)$ is asymptotically stable while the DEE point does not exist. The solution obtained by the proposed NSFD scheme (3.5) confirms the stability of the DFE point, see Figure 4.(a). For comparison, we also plot the numerical solutions obtained by the ordinary differential equations (ODE) solvers in MATLAB, namely ode23s (the $2^{nd}$ order modified Rosenbrock), ode23t (the trapezoidal rule using a "free" interpolant) and ode45 (the Runge-Kutta Dormand-Prince method of order (4, 5)) with their default tolerance. The $y$-axis (susceptible and infective sub-populations) in this figure is plotted in the logarithmic scale. Notice that the logarithmic function $\log (y)$ is defined only for $y > 0$. Hence, if the solution (susceptible or infective sub-population) is negative then the figure cannot be plotted (it is ignored by MATLAB). Figure 4.(a) shows that the solution given by the proposed NSFD scheme is convergent to the DFE point and is positive for all $t$. We observe in Figure 4.(b-d) that solutions obtained by all three MATLAB ODE solvers go to the DFE point but the infective subpopulation $I(t)$ is discontinuous at some points, showing that all three MATLAB ODE solvers produce unrealistic negative solutions. Although the tolerances of MATLAB ODE solvers are reduced such that AbsTol $= 10^{-12}$ and RelTol $= 10^{-12}$, we still observe that the three MATLAB ODE solvers produce unrealistic negative solutions.

5.   Conclusions

In this article we have proposed and analyzed a dynamically-consistent discrete SIR epidemic model with saturated incidence rate and vaccination. The discrete model is constructed using the NSFD method. In contrast to the discrete NSFD model available in literature, we constructed the discrete NSFD model by applying a nonlocal approximation for the nonlinear terms and choose a suitable denominator function. The denominator function is derived in such away that the NSFD scheme satisfies the exact boundedness of solutions. The existence conditions of all equilibrium points of the proposed discrete NSFD model are derived and the sufficient conditions for which the equilibrium points are locally and globally stable have also been established. We have shown that the proposed discrete NSFD model has exactly the same dynamics as those of continuous model irrespective of the time-step size. These qualitative properties have been confirmed by our numerical simulations. Furthermore, our numerical simulations showed that the discrete model obtained by the explicit Euler method and the NSFD scheme proposed in other literatur is dynamically-consistent with the related continuous model only when the time-step size is relatively small. The three MATLAB ODE solvers (ode23s, ode23t, ode45) may fail to produce realistic solutions, i.e., the solutions may be negative.

Acknowledgments

This research is funded by FMIPA-UB via PNBP-University of Brawijaya according to DIPA-UB No. DIPA-023.17.2.677512/2020, under contract No. 12/UN10.F09/PN/2020.

Conflict of interest

All authors declare no conflict of interest.