# Dynamical analysis of a generalized hepatitis B epidemic model and
its dynamically consistent discrete model

Manh Tuan Hoang<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of
Mathematics, FPT University, Hoa Lac Hi-Tech Park, Km29 Thang Long Blvd,
Hanoi, Viet Nam

###### Abstract

The aim of this work is to study qualitative dynamical properties of a
generalized hepatitis B epidemic model and its dynamically consistent discrete
model. Positivity, boundedness, the basic reproduction number and asymptotic
stability properties of the model are analyzed rigorously. By the Lyapunov
stability theory and the Poincare-Bendixson theorem in combination with the
Bendixson-Dulac criterion, we show that a disease-free equilibrium point is
globally asymptotically stable if the basic reproduction number
$\mathcal{R}_{0}\leq 1$ and a disease-endemic equilibrium point is globally
asymptotically stable whenever $\mathcal{R}_{0}>1$. Next, we apply the
Mickens’ methodology to propose a dynamically consistent nonstandard finite
difference (NSFD) scheme for the continuous model. By rigorously mathematical
analyses, it is proved that the constructed NSFD scheme preserves essential
mathematical features of the continuous model for all finite step sizes.
Finally, numerical experiments are conducted to illustrate the theoretical
findings and to demonstrate advantages of the NSFD scheme over standard ones.
The obtained results in this work not only improve but also generalize some
existing recognized works.

###### keywords:

HBV , Dynamical analysis , Asymptotic stability , NSFD schemes , Dynamic
consistency

###### MSC:

[2020] 34C60, 34D05, 37C75, 37N99

††journal: Mathematical Methods in the Applied Sciences

##  1 Introduction

Mathematical modeling and analysis of infectious diseases has been a central
research field in infectious disease epidemiology with many useful
applications not only in theory but also in practice [3, 7, 8, 9, 27, 28, 35,
41, 44, 55]. The well-known SIR model, proposed by Kermack and McKendrick in
1927 [35], can be considered as one of the first epidemic models and is
usually used to introduce epidemic modeling. Over the past of several decades,
for modeling purpose, a great number of mathematical models have been studied
to predict and discover transmission mechanisms of infectious diseases. As an
important consequence, strategies and measures for controlling and preventing
infectious diseases can be provided.

It is well-known that hepatitis B is a dangerous infectious disease that
attacks the liver and can cause acute and chronic diseases. Nowadays,
hepatitis B has been a major global health problem. Hence, effective
strategies and measures for preventing and controlling hepatitis B virus (HBV)
are always needed. For this reason, a large number of mathematical models have
been proposed and studied by many mathematicians and epidemiologists (see, for
instance, [2, 11, 12, 14, 21, 22, 24, 36, 37, 38, 43, 58]). These mathematical
models can help us predict and discover transmission mechanisms and
characteristics of the HBV. Based on this, strategies and measures for
projecting public health and for preventing the hepatitis B can be suggested.

In [36], Khan et al. proposed a hepatitis B epidemic model based on HBV
transmission characteristics. The model is described by a system of nonlinear
differential equations

$\begin{split}&\dfrac{dS(t)}{dt}=\Lambda-\dfrac{\alpha S(t)I(t)}{1+\gamma
I(t)}-\big{(}\mu_{0}+\nu\big{)}S(t),\\\ &\dfrac{dI(t)}{dt}=\dfrac{\alpha
S(t)I(t)}{1+\gamma I(t)}-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I(t),\\\
&\dfrac{dR(t)}{dt}=\beta I(t)+\nu S(t)-\mu_{0}R(t),\end{split}$ (1)

subject to initial data

$S(0)>0,\quad I(0)\geq 0,\quad R(0)>0.$

In this model:

  * 1.

the entire population, denoted by $N$, is divided into $3$ classes:
susceptible ($S$), infected ($I$) and recovered or removed ($R$) classes;

  * 2.

$\Lambda$ is the birth rate;

  * 3.

$\lambda$ is the transmission rate;

  * 4.

$\mu_{0}$ and $\mu_{1}$ denote the natural and disease induced death rates,
respectively;

  * 5.

$\beta$ is the recovery rate;

  * 6.

$\nu$ and $\gamma$ stand for the vaccination and saturation rates,
respectively;

We refer the readers to [36] for more details of the model (1) and its
qualitative dynamical properties. In [32], A study on the global stability and
numerical solutions of the model (1) was performed. Besides, Hoang and
Egbelowo in [29] proposed and analyzed a generalized version of the model (1)
using the Caputo fractional derivative.

In the model (1), let us denote by

$f(S,I)=\dfrac{\alpha S(t)I(t)}{1+\gamma I(t)}.$ (2)

It is easy to verify that this function possesses the following properties:
(H1): $f:\mathbb{R}_{+}^{2}\to\mathbb{R}_{+}$ is a differentiable function,
$f(S,0)=f(0,I)=0$ for all $S,I\geq 0$ and $f(S,I)>0$ for all $S,I>0$;
(H2): there exists $\eta>0$ such that $f(S,I)\leq\eta S$ for all $S,I\geq 0$;
(H3): $\dfrac{\partial f(S,I)}{\partial S}>0$ and bounded for all $S\geq 0$
and $I>0$;
(H4): $\dfrac{\partial f(S,I)}{\partial I}\geq 0$ for all $S,I\geq 0$;
(H5): $I\dfrac{\partial f(S,I)}{\partial I}-f(S,I)\leq 0$ for all $S,I\geq 0$.
The conditions (H1)-(H5) are biologically motivated and can be found in
previous works [34, 46, 60]. In a recent work [34], Karaji and Nyamoradi
considered the model (1) in the context of the Caputo fractional-order
derivative and general incidence functions. It should be emphasized that
functions $f$ satisfying the conditions (H1)-(H5) include many famous
incidence functions, for examples

$f(S,I)=\beta SI,\quad f(S,I)=\dfrac{\beta SI}{1+bI},\quad f(S,I)=\dfrac{\beta
SI}{1+aS+bI},\quad f(S,I)=\dfrac{\beta SI}{1+aS+bI+cSI}.$

Motivated and inspired by the recognized works [34, 46, 60] as well as the
importance of mathematical models of the HBV, in this work we consider a
generalized version of the model (1) by replacing the function $f$ defined by
(2) by general ones $f$ satisfying the conditions (H1)-(H5). More precisely,
we investigate the following model

$\begin{split}&\dfrac{dS(t)}{dt}=\Lambda-f(S(t),I(t))-\big{(}\mu_{0}+\nu\big{)}S(t),\\\
&\dfrac{dI(t)}{dt}=f(S(t),I(t))-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I(t),\\\
&\dfrac{dR(t)}{dt}=\beta I(t)+\nu S(t)-\mu_{0}R(t),\end{split}$ (3)

where $f(S,I)$ is any incidence function satisfying the conditions (H1)-(H5).
Our main objective is to study qualitative dynamical properties and reliable
numerical solutions for the model (3).

In the first part of this work, we analyze positivity, boundedness, the basic
reproduction number, equilibria (disease-free and disease-endemic equilibrium
points) and asymptotic stability properties of the model (3). By the Lyapunov
stability theory and the Poincare-Bendixson theorem in combination with the
Bendixson-Dulac criterion, we prove that the disease-free equilibrium (DFE)
point is globally asymptotically stable if the basic reproduction number
$\mathcal{R}_{0}\leq 1$; and the disease-endemic equilibrium (DEE) point
exists and is globally asymptotically stable whenever $\mathcal{R}_{0}>1$.
Consequently, qualitative dynamical properties of the model (3) are fully
determined.

In the second part, we apply the Mickens’ methodology [48, 49, 50, 51, 52] to
propose a dynamically consistent nonstandard finite difference (NSFD) scheme
for the model (3). By rigorously mathematical analyses, we show that the
proposed NSFD scheme preserves essential mathematical features of the
continuous model for all finite step sizes. In other words, the NSFD scheme
behaves similarly to the continuous model regardless of chosen step sizes. It
is worth noting that one of the prominent advantages of NSFD schemes is that
they are able to correctly maintain important mathematical features of
solutions of differential equations models (positivity, boundedness,
monotonicity, stability, periodicity, physical properties, etc.) for any
finite step size. Hence, NSFD schemes are effective to simulate dynamics of
dynamical differential models over long time periods. Because of this,
nowadays NSFD schemes have been recognized as one of the effective approaches
to solve differential equation equations arising in theory and practice [1, 5,
10, 13, 23, 25, 26, 39, 53, 56, 57, 62, 63]. In recent works [15, 16, 17, 18,
19, 20, 31, 32], we have successfully developed the Mickens’s methodology to
construct NSFD schemes for some mathematical models arising in biology and
epidemiology.

In the third part, we conduct numerical experiments to illustrate the
theoretical findings and to demonstrate advantages of the NSFD scheme over
standard ones. The numerical results show that the Euler and second-order
Runge-Kutta (RK2) schemes can generate numerical approximations which destroy
not only the positivity but also the asymptotic stability properties of the
continuous model for some given step sizes; meanwhile, the NSFD scheme
preserve these properties for the same step sizes.

It is worth noting that the mathematical analyses in [36] failed to conclude
the global asymptotic stability (GAS) of the DEE point of the model (1).
However, we obtain the GAS of the model (1) thanks to the global stability
analysis for the model (3) (Section 2). This provides an important improvement
for the results in [36]. In a recent work [59], Suryanto and Darti constructed
an NSFD scheme for the model (1) but the convergence and error bounds for this
NSFD scheme were not established. Hence, our constructed NSFD scheme not only
generalizes but also provides the convergence and error estimates for the
Suryanto and Darti’s scheme.

The plan of this work is as follows:
Dynamical properties of the model (3) are studied in Section 2. The
dynamically consistent NSFD scheme is formulated and analyzed in Section 3. In
Section 4, we report some numerical experiments to illustrate the theoretical
results. Some discussions and remarks are provided in the last section.

##  2 Dynamical analysis

In this section, dynamical qualitative properties of the model (3) will be
investigated. We first establish the positivity and boundedness of the model.

######  Proposition 1 (Positivity and boundedness).

The set $\mathbb{R}_{+}^{3}$ is a positively invariant set of the model (3),
i.e., $S(t),I(t),R(t)\geq 0$ whenever $S(0),I(0),R(0)\geq 0$. Furthermore, we
have

$\begin{split}&\limsup_{t\to\infty}S(t)\leq\dfrac{\Lambda}{\mu_{0}+\nu},\quad\dfrac{\Lambda}{p}\leq\limsup_{t\to\infty}\big{(}S(t)+I(t)\big{)}\leq\dfrac{\Lambda}{q},\quad\dfrac{r_{1}\Lambda}{p\mu_{0}}\leq\limsup_{t\to\infty}R(t)\leq\dfrac{r_{2}\Lambda}{q\mu_{0}},\\\
&p:=\max\\{\mu_{0}+\nu,\,\,\mu_{0}+\mu_{1}+\beta\\},\quad
q:=\min\\{\mu_{0}+\nu,\,\,\mu_{0}+\mu_{1}+\beta\\},\quad
r_{2}:=\max\\{\beta,\,\,\nu\\},\quad
r_{1}:=\min\\{\beta,\,\,\nu\\}.\end{split}$ (4)

###### Proof.

First, the system (3) implies that

$\dfrac{dS}{dt}\bigg{|}_{S=0}=\Lambda>0,\quad\dfrac{dI}{dt}\bigg{|}_{I=0}=0,\quad\dfrac{dR}{dt}\bigg{|}_{R=0}=\beta
I+\nu S\geq 0.$

Thanks to [33, Lemma 10], the positivity of the solutions of the model (3) is
obtained. Hence, the set $\mathbb{R}_{3}^{+}$ is a positively invariant set of
(3).

Next, it follows from the first equation of (3) that

$\dfrac{dS}{dt}\leq\Lambda-(\mu_{0}+\nu)S.$

Using the basic comparison theorem for ODEs [47], we obtain

$S(t)\leq\bigg{(}S(0)-\dfrac{\Lambda}{\mu_{0}+\nu}\bigg{)}e^{-(\mu_{0}+\nu)t}+\dfrac{\Lambda}{\mu_{0}+\nu},$

which implies the first estimate of (4).

Similarly, it follows from the first and second equation of (3) that

$\Lambda-p(S+I)\leq\dfrac{d(S+I)}{dt}=\Lambda-(\mu_{0}+\nu)S-(\mu_{0}+\mu_{1}+\beta)I\leq\Lambda-q(S+I).$

Hence, the second estimate of (4) is obtained.

Lastly, using the third equation of (3) we have

$\dfrac{r_{1}\Lambda}{p}-\mu_{0}R\leq
r_{1}(S+I)-\mu_{0}R\leq\dfrac{dR}{dt}\leq
r_{2}(S+I)-\mu_{0}R\leq\dfrac{r_{2}\Lambda}{q}-\mu_{0}R$

for $t$ large enough. Consequently, the third estimate of (4) is proved. This
proof is complete. ∎

######  Proposition 2 (Equilibria and the basic reproduction number).

  1. (i)

The HBV model (3) always possesses a disease-free equilibrium (DFE) point
given by

$E^{0}=(S^{0},I^{0},R^{0})=\bigg{(}\dfrac{\Lambda}{\mu_{0}+\nu},\,\,0\,\,\dfrac{\nu\Lambda}{\mu_{0}(\mu_{0}+\nu)}\bigg{)}$
(5)

for all the values of the parameters.

  2. (ii)

The basic reproduction number of the model (3) can be computed as

$\mathcal{R}_{0}=\dfrac{1}{\mu_{0}+\mu_{1}+\beta}\dfrac{\partial
f\bigg{(}\dfrac{\Lambda}{\mu_{0}+\nu},\,\,0\bigg{)}}{\partial I}.$ (6)

  3. (iii)

The model (3) has a unique disease-endemic equilibrium (DEE) point
$E^{*}=(S^{*},I^{*},R^{*})$ if and only if $\mathcal{R}_{0}>1$. Moreover, if
$E^{*}$ exists it is given by

$S^{*}\dfrac{\Lambda-(\mu_{0}+\mu_{1}+\beta)I^{*}}{\mu_{0}+\nu},\quad
R^{*}=\dfrac{\beta I^{*}+\nu S^{*}}{\mu_{0}},$ (7)

where $I^{*}$ is a unique positive solution of the equation

${f\bigg{(}\dfrac{\Lambda-(\mu_{0}+\mu_{1}+\beta)I}{\mu_{0}+\nu},\,\,I\bigg{)}}-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I=0.$

###### Proof.

Proof of Part (i). To determine equilibria of the model (3), we need to solve
the nonlinear system of algebra equations

$\Lambda-f(S,I)-\big{(}\mu_{0}+\nu\big{)}S=0,\quad
f(S,I)-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I=0,\quad\beta I+\nu S-\mu_{0}R=0.$
(8)

It follows from the first and second equations of the system (8) that

$S=\dfrac{\Lambda-(\mu_{0}+\mu_{1}+\beta)I}{\mu_{0}+\nu}.$

Combining this with the second equation of (8) we obtain

$f\bigg{(}\dfrac{\Lambda-(\mu_{0}+\mu_{1}+\beta)I}{\mu_{0}+\nu},I\bigg{)}-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I=0.$
(9)

Since $f(S,0)=0$ for all $S,I\geq 0$, the equation (9) always has a trivial
solution $I=0$. Consequently, the system (8) always possesses a solution
$\big{(}\Lambda/(\mu_{0}+\nu),\,\,0,\,\,{\nu\Lambda}/{[\mu_{0}(\mu_{0}+\nu)]}\big{)}$.
This solution corresponds to a unique DFE point.
Proof of Part (ii). We apply the method proposed in [61] to compute the basic
reproduction number of the model (3). Since the first two equations of (3) do
not depend on $R$, it is sufficient to consider the following sub-model

$\begin{split}&\dfrac{dS(t)}{dt}=\Lambda-f(S(t),I(t))-\big{(}\mu_{0}+\nu\big{)}S(t),\\\
&\dfrac{dI(t)}{dt}=f(S(t),I(t))-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I(t).\end{split}$
(10)

We reorder the variables in (10) as $(I,S)$. Then, the DFE point is
transformed to $x_{0}=(I_{0},S_{0})$. The model (10) can be written in the
matrix form

$\dfrac{dx}{dt}=\mathcal{F}(x)-\mathcal{V}(x),$

where

$\mathcal{F}(x)=\begin{pmatrix}f(S,I)\\\
\Lambda\end{pmatrix},\quad\mathcal{V}(x)=\begin{pmatrix}(\mu_{0}+\mu_{1}+\beta)I\\\
f(S,I)+(\mu_{0}+\nu)S\end{pmatrix}$

Consequently,

$D\mathcal{F}(x_{0})=\begin{pmatrix}\dfrac{\partial f}{\partial
I}(S_{0},0)&0\\\ &\\\ 0&0\end{pmatrix},\quad
D\mathcal{V}(x)=\begin{pmatrix}(\mu_{0}+\mu_{1}+\beta)&0\\\ &\\\
\dfrac{\partial f}{\partial I}(S_{0},0)&\dfrac{\partial f}{\partial
S}(S_{0},0)+(\mu_{0}+\nu).\end{pmatrix}$

Then, the basic reproduction number can be computed as

$\mathcal{R}_{0}=\rho(FV^{-1})=\dfrac{1}{\mu_{0}+\mu_{1}+\beta}\dfrac{\partial
f\bigg{(}\dfrac{\Lambda}{\mu_{0}+\nu},\,\,0\bigg{)}}{\partial I}.$

Proof of Part (iii). Return to the equation (9). We will prove that it has
only a unique positive solution $I^{*}>0$. Indeed, consider the equation

$F(I):=\dfrac{f\bigg{(}\dfrac{\Lambda-(\mu_{0}+\mu_{1}+\beta)I}{\mu_{0}+\nu},\,\,I\bigg{)}}{I}-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}=0.$
(11)

From the condition (H5), it is easy to verify that $F^{\prime}(I)<0$ for all
$I>0$. On the other hand,

$\begin{split}&\lim_{I\to 0}F(I)=\dfrac{\partial
f\bigg{(}\dfrac{\Lambda}{\mu_{0}+\nu},\,\,0\bigg{)}}{\partial
I}-({\mu_{0}+\mu_{1}+\beta})>0,\\\
&F\bigg{(}\dfrac{\Lambda}{\mu_{0}+\mu_{1}+\beta}\bigg{)}=-({\mu_{0}+\mu_{1}+\beta})<0.\end{split}$

Hence, the equation (11) has a positive solution $I^{*}$ and we obtain
$(S^{*},I^{*},R^{*})$ as in (7). In other words, the existence and uniqueness
of the DEE point are proved. ∎

######  Lemma 1 (Local asymptotic stability).

  1. (i)

The DFE point $E^{0}$ is locally asymptotically stable if $\mathcal{R}_{0}<1$
and unstable if $\mathcal{R}_{0}>1$.

  2. (ii)

The DEE point $E^{*}$ is locally asymptotically stable if $\mathcal{R}_{0}>1$
and unstable if $\mathcal{R}_{0}<1$.

###### Proof.

Proof of Part (i). The Jacobian matrix of the system (3) at $E^{0}$ is given
by

$J(E^{0})=\begin{pmatrix}-(\mu_{0}+\nu)&-\dfrac{\partial f(E^{0})}{\partial
I}&0\\\ &&\\\ 0&\dfrac{\partial f(E^{0})}{\partial
I}-(\mu_{0}+\mu_{1}+\beta)&0\\\ &&\\\ \nu&\beta&-\mu_{0}\end{pmatrix}.$

Hence, three eigenvalues of $J(E^{0})$ are

$\lambda_{1}=-(\mu_{0}+\nu),\quad\lambda_{2}=\dfrac{\partial
f(E^{0})}{\partial I}-(\mu_{0}+\mu_{1}+\beta),\quad\lambda_{3}=-\mu_{0}.$

It is clear that $\lambda_{1}<0$ and $\lambda_{3}<0$. Moreover, if
$\mathcal{R}_{0}<1$ then $\lambda_{3}<0$. By the Lyapunov direct method [3,
44], the DFE point $E^{0}$ is locally asymptotically stable. Conversely, if
$\mathcal{R}_{0}>1$ then $\lambda_{3}>0$ and hence, $E^{0}$ is unstable. The
proof is complete.
Proof of Part (ii). Note that $E^{*}$ exists if and only if
$\mathcal{R}_{0}>1$. The Jacobian matrix of the system (3) at $E^{*}$ is given
by

$J(E^{*})=\begin{pmatrix}-\dfrac{\partial f(E^{*})}{\partial
S}-(\mu_{0}+\nu)&-\dfrac{\partial f(E^{*})}{\partial I}&0\\\ &&\\\
\dfrac{\partial f(E^{*})}{\partial S}&\dfrac{\partial f(E^{*})}{\partial
I}-(\mu_{0}+\mu_{1}+\beta)&0\\\ &&\\\ \nu&\beta&-\mu_{0}\end{pmatrix}.$

Hence, one of three eigenvalues of $J(E^{*})$ is $\lambda_{1}=-\mu_{0}<0$ and
two remaining eigenvalues are the ones of the sub-matrix

$\widehat{J}(E^{*})=\begin{pmatrix}-\dfrac{\partial f(E^{*})}{\partial
S}-(\mu_{0}+\nu)&-\dfrac{\partial f(E^{*})}{\partial I}\\\ &&\\\
\dfrac{\partial f(E^{*})}{\partial S}&\dfrac{\partial f(E^{*})}{\partial
I}-(\mu_{0}+\mu_{1}+\beta)\end{pmatrix}.$

Because $E^{*}$ is the positive equilibrium point we have

$\dfrac{f(S^{*},I^{*})}{I^{*}}=(\mu_{0}+\mu_{1}+\beta).$

Hence, it follows from the hypothesis (H5) that

$\dfrac{\partial f(S^{*},I^{*})}{\partial
I}<\dfrac{f(S^{*},I^{*})}{I^{*}}=(\mu_{0}+\mu_{1}+\beta).$ (12)

From this estimate, we obtain

$Trace(\widehat{J})<0,\quad\det(\widehat{J})>0.$

Thanks to the Lyapunov stability theorem [3, 44], the local asymptotic
stability of $E^{*}$ is proved. This proof is completed. ∎

To end this section, we investigate the GAS of the model (3). As an
consequence of Proposition 1, it is sufficient to consider the GAS of the
model (3) on its feasible set defined by

$\Omega=\bigg{\\{}(S,I,R)\Big{|}S,I,R\geq 0,\,\,0\leq
S\leq\dfrac{\Lambda}{\mu_{0}+\nu},\,\,\dfrac{\Lambda}{p}\leq
S+I\leq\dfrac{\Lambda}{q},\,\,\dfrac{r_{1}\Lambda}{p\mu_{0}}\leq
R\leq\dfrac{r_{2}\Lambda}{q\mu_{0}}\bigg{\\}},$

where $p,q,r_{1}$ and $r_{2}$ are given by (4). Moreover, since the first two
equations of the model (3) do not depend on $R$, we only need to consider the
GAS of the sub-model

$\begin{split}&\dfrac{dS}{dt}=\Lambda-f(S,I)-\big{(}\mu_{0}+\nu\big{)}S,\\\
&\dfrac{dI}{dt}=f(S,I)-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I\end{split}$ (13)

on its feasible set

$\Omega^{*}=\bigg{\\{}(S,I)\big{|}S,I\geq 0,\,\,0\leq
S\leq\dfrac{\Lambda}{\mu_{0}+\nu},\,\,\dfrac{\Lambda}{p}\leq
S+I\leq\dfrac{\Lambda}{q}\bigg{\\}}.$

On the set $\Omega^{*}$, the model (13) always has a DFE point
$\widehat{E}^{0}=(S^{0},0)$ for all the values of the parameters and a DEE
point $\widehat{E}^{*}=(S^{*},I^{*})$ exists if and only if
$\mathcal{R}_{0}>1$.

######  Lemma 2.

Let $f(S,I)$ be a function satisfying the conditions (H1)-(H5). We define

$V(S):=S-S_{0}-\int_{S_{0}}^{S}\dfrac{f(S_{0},I)}{f(t,I)}dt,\quad S\geq 0$
(14)

and

$g(I)=\dfrac{f(S_{0},I)}{I}-(\mu_{0}+\mu_{1}+\nu),\quad I>0.$ (15)

Then,

  1. (i)

$V(S)\geq 0$ for all $S\geq 0$ and $V(S)=0$ if and only if $S=S_{0}$.

  2. (ii)

If $\mathcal{R}_{0}\leq 1$, then $g(I)<0$ for all $I>0$.

###### Proof.

Proof of Part (i). First, it is clear that $V(S_{0})=0$. The derivative of the
function $V$ is given by

$V^{\prime}(S)=1-\dfrac{f(S_{0},I)}{f(S,I)}.$

From the condition (H3) we deduce that

$\begin{split}&V^{\prime}(S)>0\quad\mbox{for all}\quad S>S_{0},\\\
&V^{\prime}(S)<0\quad\mbox{for all}\quad S<S_{0},\\\
&V^{\prime}(S)=0\quad\mbox{if and only if}\quad S=S_{0}.\end{split}$

Hence, $V(S)>V(S_{0})=0$ for all $S\neq S_{0}$. The proof of this part is
complete.
Proof of Part (ii). It follows from the condition (H5) that

$g^{\prime}(I)=\dfrac{\partial f(S_{0},I)}{\partial I}<0.$

On the other hand,

$\lim_{I\to 0}g(I)=\dfrac{\partial
f\bigg{(}\dfrac{\Lambda}{\mu_{0}+\nu},\,\,0\bigg{)}}{\partial
I}-(\mu_{0}+\mu_{1}+\nu)=(\mu_{0}+\mu_{1}+\nu)(\mathcal{R}_{0}-1)\leq 0.$

Therefore, we conclude that $g(I)<(\mu_{0}+\mu_{1}+\nu)(\mathcal{R}_{0}-1)\leq
0$. This completes the proof. ∎

######  Theorem 1.

The DFE point $\widehat{E^{0}}$ of the model (13) is globally asymptotically
stable if $\mathcal{R}_{0}\leq 1$.

###### Proof.

Consider a Lyapunov function given by

$L(S,I)=V(S)+I,$ (16)

where $V(S)$ is defined by (14). The derivative of the function $L$ along
solutions of (13) is given by

$\begin{split}\dfrac{dL}{dt}=\dfrac{dL}{dS}\dfrac{dS}{dt}+\dfrac{dL}{dI}\dfrac{dI}{dt}&=\bigg{(}1-\dfrac{f(S_{0},I)}{f(S,I)}\bigg{)}\Big{[}\Lambda-f(S,I)-\big{(}\mu_{0}+\nu\big{)}S\Big{]}+\Big{[}f(S,I)-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I\Big{]}\\\
&=\bigg{(}1-\dfrac{f(S_{0},I)}{f(S,I)}\bigg{)}\Big{[}\Lambda-\big{(}\mu_{0}+\nu\big{)}S\Big{]}+\Big{[}f(S_{0},I)-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I\Big{]}\\\
&=\bigg{(}1-\dfrac{f(S_{0},I)}{f(S,I)}\bigg{)}\big{(}\mu_{0}+\nu\big{)}\big{(}S_{0}-S\big{)}+Ig(I),\end{split}$

where the function $g(I)$ is defined by (15).

Combining the properties of the function $f(S,I)$ and Lemma 2, we deduce that
the function $L$ given by (16) satisfies the Lyapunov direct method [40, 42].
So, the GAS of $\widehat{E^{0}}$ is obtained. The proof is complete. ∎

######  Theorem 2.

Suppose that $\mathcal{R}_{0}>1$. Then, the DEE point $\widehat{E^{*}}$ of the
model (13) is globally asymptotically stable whenever $I(0)>0$.

###### Proof.

Note that if $I(0)=0$ then $\big{(}S(t),I(t)\big{)}$ will converge to the DFE
point $\widehat{E^{0}}$ as $t\to\infty$. We first use the Dulac-Bendixson
criterion [44, Theorem 3.6] to show that the system (13) has no periodic
orbits or graphics in $\mathbb{R}_{2}^{+}$. Indeed, let us denote by $F(S,I)$
and $G(S,I)$ be the right-hand side functions of (13). Consider a function

$D(S,I)=\dfrac{1}{I}.$

We have

$\begin{split}\dfrac{\partial(DF)}{\partial S}+\dfrac{\partial(DG)}{\partial
I}&=\dfrac{\partial}{\partial
S}\bigg{[}\dfrac{\Lambda}{I}-\dfrac{f(S,I)}{I}-\dfrac{(\mu_{0}+\nu)S}{I}\bigg{]}+\dfrac{\partial}{\partial
I}\bigg{[}\dfrac{f(S,I)}{I}-(\mu_{0}+\mu_{1}+\beta)\bigg{]}\\\
&=-\dfrac{1}{I}\dfrac{\partial
f(S,I)}{S}-\dfrac{\mu_{0}+\nu}{I}+\dfrac{\partial}{\partial
I}\bigg{(}\dfrac{f(S,I)}{I}\bigg{)}.\end{split}$

It follows from the conditions (H3) and (H5) that

$\dfrac{\partial(DF)}{\partial S}+\dfrac{\partial(DG)}{\partial I}<0.$

Hence, the system (13) has no periodic orbits or graphics in the open first
quadrant.

Since all solutions of the model (13) are bounded and the DFE point
$\widehat{E^{0}}$ is unstable saddle if $\mathcal{R}_{0}>1$, using Poincaré-
Bendixson theorem [44, Theorem 3.5] we obtain the GAS of the DEE point
$\widehat{E^{*}}$. The proof is complete. ∎

From Theorems 1 and 2, we obtain the GAS of the model (3) as follows.

######  Theorem 3.

  1. (i)

The DFE point ${E^{0}}$ of the model (3) is not only locally asymptotically
stable but also globally asymptotically stable if $\mathcal{R}_{0}\leq 1$.

  2. (ii)

The DEE point ${E^{*}}$ of the model (3) is not only locally asymptotically
but also globally asymptotically stable provided that $\mathcal{R}_{0}>1$ and
$I(0)>0$.

##  3 Construction of dynamically consistent NSFD scheme

In this section, we construct and analyze a dynamically consistent NSFD scheme
for the model (3). For this purpose, consider the model (3) on a finite time
interval $[0,T]$ and discretize this interval by

$0=t_{0}<t_{1}<\ldots<t_{N-1}<t_{N}=T,$

where $t_{n+1}-t_{n}=\Delta t$ for $0\leq n\leq N-1$. Let us denote by
$(S_{n},I_{n},R_{n})$ the intended approximation for
$\big{(}S(t_{n}),I(t_{n}),R(t_{n})\big{)}$ for $n=1,2,\ldots,N$. Applying the
Mickens’ methodology [48, 49, 50, 51, 52] we approximate the first derivatives
in the model (3) by

$\dfrac{dS}{dt}\bigg{|}_{t=t_{n}}\approx\dfrac{S_{n+1}-S_{n}}{\varphi(\Delta
t)},\qquad\dfrac{dI}{dt}\bigg{|}_{t=t_{n}}\approx\dfrac{I_{n+1}-I_{n}}{\varphi(\Delta
t)},\qquad\dfrac{dR}{dt}\bigg{|}_{t=t_{n}}\approx\dfrac{R_{n+1}-R_{n}}{\varphi(\Delta
t)},$ (17)

where $\varphi$, called the denominator function, is a positive function
satisfying $\varphi(\Delta t)=\Delta t+\mathcal{O}(\Delta t^{2})$ as $\Delta
t\to 0$. Next, we discretize the right-hand side function of the model (3) as
follows. Let us denote

$f^{*}(S,I)=\begin{cases}&\dfrac{f(S,I)}{S},\quad S>0,\\\ &0,\quad
S=0.\end{cases}$ (18)

Then, $f(S,I)=Sf^{*}(S,I)$. Moreover, it follows from the conditions (H1) and
(H5) that $f^{*}(S,I)$ is differentiable for $S,I\geq 0$. In particular,

$\begin{split}&f^{*}(S^{0},I^{0})=f^{*}(S^{*},I^{*})=0,\\\ &\dfrac{\partial
f^{*}}{\partial S}(S_{0},0)=0,\quad\dfrac{\partial f^{*}}{\partial
I}(S_{0},0)=\dfrac{1}{S_{0}}\dfrac{\partial f}{\partial I}(S_{0},0),\\\
&\dfrac{\partial f^{*}}{\partial
S}(S^{*},I^{*})=\dfrac{1}{S^{*}}\dfrac{\partial f}{\partial
S}(S^{*},I^{*}),\quad\dfrac{\partial f^{*}}{\partial
I}(S^{*},I^{*})=\dfrac{1}{S^{*}}\dfrac{\partial f}{\partial
I}(S^{*},I^{*}).\end{split}$ (19)

We now discretize the right-hand side function of (3) by

$\begin{split}&\Lambda-f(S(t_{n}),I(t_{n}))-\big{(}\mu_{0}+\nu\big{)}S(t_{n})\approx\Lambda-
S_{n+1}f^{*}(S_{n},I_{n})-\big{(}\mu_{0}+\nu\big{)}S_{n+1}\\\
&f(S(t_{n}),I(t_{n}))-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I(t_{n})\approx
S_{n+1}f^{*}(S_{n},I_{n})-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I_{n+1},\\\
&\beta I(t_{n})+\nu S(t_{n})-\mu_{0}R(t_{n})\approx\beta I_{n+1}+\nu
S_{n+1}-\mu_{0}R_{n+1}.\end{split}$ (20)

The systems (17) and (20) lead to the the following NSFD model for the model
(3)

$\begin{split}&\dfrac{S_{n+1}-S_{n}}{\varphi(\Delta t)}=\Lambda-
S_{n+1}f^{*}(S_{n},I_{n})-\big{(}\mu_{0}+\nu\big{)}S_{n+1},\\\
&\dfrac{I_{n+1}-I_{n}}{\varphi(\Delta
t)}=S_{n+1}f^{*}(S_{n},I_{n})-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I_{n+1},\\\
&\dfrac{R_{n+1}-R_{n}}{\varphi(\Delta t)}=\beta I_{n+1}+\nu
S_{n+1}-\mu_{0}R_{n+1}.\end{split}$ (21)

######  Theorem 4.

The set $\mathbb{R}^{3}_{+}$ is a positively invariant set of the NSFD model
(3), that is, $S_{n},I_{n},R_{n}\geq 0$ whenever $S_{0},I_{0},R_{0}\geq 0$.
Moreover, we have

$\begin{split}&\limsup_{n\to\infty}S_{n}\leq\dfrac{\Lambda}{\mu_{0}+\nu},\quad\quad\dfrac{\Lambda}{p}\leq\limsup_{n\to\infty}\big{(}S_{n}+I_{n}\big{)}\leq\dfrac{\Lambda}{q},\quad\dfrac{r_{1}\Lambda}{p\mu_{0}}\leq\limsup_{n\to\infty}R_{n}\leq\dfrac{r_{2}\Lambda}{q\mu_{0}},\end{split}$
(22)

where $p,q,r_{1}$ and $r_{2}$ are given by (4).

###### Proof.

First, we rewrite the system (21) in the explicit form

$\begin{split}&S_{n+1}=\dfrac{S_{n}+\varphi\Lambda}{1+\varphi
f^{*}(S_{n},I_{n})+\varphi(\mu_{0}+\nu)},\\\ &I_{n+1}=\dfrac{I_{n}+\varphi
S_{n+1}f^{*}(S_{n},I_{n})}{1+\varphi(\mu_{0}+\mu_{1}+\beta)},\\\
&R_{n+1}=\dfrac{R_{n}+\varphi\beta I_{n+1}+\varphi\nu
S_{n+1}}{1+\varphi\mu_{0}}.\end{split}$ (23)

Hence, $S_{n+1},I_{n+1},R_{n+1}\geq 0$ if $S_{n},I_{n},R_{n}\geq 0$.

It follows from the first equation of (21) that

$\dfrac{S_{n+1}-S_{n}}{\varphi}\leq\Lambda-(\mu_{0}+\nu)S_{n+1},$

which implies that

$\begin{split}S_{n+1}&\leq\dfrac{S_{n}}{1+\varphi(\mu_{0}+\nu)}+\dfrac{\varphi\Lambda}{1+\varphi(\mu_{0}+\nu)}\\\
&\leq\dfrac{1}{1+\varphi(\mu_{0}+\nu)}\bigg{[}\dfrac{S_{n-1}}{1+\varphi(\mu_{0}+\nu)}+\dfrac{\varphi\Lambda}{1+\varphi(\mu_{0}+\nu)}\bigg{]}+\dfrac{\varphi\Lambda}{1+\varphi(\mu_{0}+\nu)}\\\
&=\bigg{[}\dfrac{1}{1+\varphi(\mu_{0}+\nu)}\bigg{]}^{2}S_{n-1}+\dfrac{\varphi\Lambda}{1+\varphi(\mu_{0}+\nu)}\bigg{[}1+\dfrac{1}{1+\varphi(\mu_{0}+\nu)}\bigg{]}\\\
&\leq\ldots\leq\bigg{[}\dfrac{1}{1+\varphi(\mu_{0}+\nu)}\bigg{]}^{n+1}S_{0}+\dfrac{\varphi\Lambda}{1+\varphi(\mu_{0}+\nu)}\sum_{j=0}^{n}\bigg{(}\dfrac{1}{1+\varphi(\mu_{0}+\nu)}\bigg{)}^{j}\\\
&=\bigg{[}\dfrac{1}{1+\varphi(\mu_{0}+\nu)}\bigg{]}^{n+1}S_{0}+\dfrac{\varphi\Lambda}{1+\varphi(\mu_{0}+\nu)}\dfrac{1-\bigg{(}\dfrac{1}{1+\varphi(\mu_{0}+\nu)}\bigg{)}^{n+1}}{1-\dfrac{1}{1+\varphi(\mu_{0}+\nu)}}.\end{split}$

Letting $n\to\infty$, the first estimate of (22) is obtained.

Similarly, from the first and second equations of (21) we have

$\Lambda-p(S_{n+1}+I_{n+1})\leq\dfrac{(S_{n+1}+I_{n+1})-(S_{n}+I_{n})}{\varphi}=\Lambda-(\mu_{0}+\nu)S_{n+1}-(\mu_{0}+\mu_{1}+\beta)I_{n+1}\leq\Lambda-q(S_{n+1}+I_{n+1}),$

which implies that

$\dfrac{\varphi\Lambda+(S_{n}+I_{n})}{1+\varphi
p}\leq(S_{n+1}+I_{n+1})\leq\dfrac{\varphi\Lambda+(S_{n}+I_{n})}{1+\varphi q}.$

From this estimate, we obtain the second estimate of (22).

Lastly, it follows from the third equation of (21) that

$r_{1}\dfrac{\Lambda}{p}-\mu_{0}R_{n+1}\leq
r_{1}(S_{n+1}+I_{n+1})-\mu_{0}R_{n+1}\leq\dfrac{R_{n+1}-R_{n}}{\varphi}\leq
r_{2}(S_{n+1}+I_{n+1})-\mu_{0}R_{n+1}\leq
r_{2}\dfrac{\Lambda}{q}-\mu_{0}R_{n+1}$

for $n$ large enough. Hence,

$\dfrac{R_{n}+\varphi r_{1}{\Lambda}/{p}}{1+\varphi\mu_{0}}\leq
R_{n+1}\leq\dfrac{R_{n}+\varphi r_{2}{\Lambda}/{q}}{1+\varphi\mu_{0}},$

which implies the last estimate of (22). The proof is completed. ∎

We now compute the basic reproduction number for the discrete model (21) by
the next generation matrix approach [4].

It is easy to verify that the model (21) always has a unique DFE point
$E^{0}=(S_{0},0,R_{0})$ for all the values of the parameters. Since the first
two equations of (21) do not include $R_{n}$, it is sufficient to consider the
following sub-model

$\begin{split}&\dfrac{S_{n+1}-S_{n}}{\varphi(\Delta t)}=\Lambda-
S_{n+1}f^{*}(S_{n},I_{n})-\big{(}\mu_{0}+\nu\big{)}S_{n+1},\\\
&\dfrac{I_{n+1}-I_{n}}{\varphi(\Delta
t)}=S_{n+1}f^{*}(S_{n},I_{n})-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I_{n+1}.\end{split}$
(24)

We reorder the variables in (24) as $(I_{n},S_{n})$. Then, the DFE point is
transformed to $(0,S_{0})$. By using (19), the corresponding Jacobian matrix
of (24) at $(0,S_{0})$ is

$J_{0}=\begin{pmatrix}\dfrac{1+\varphi\dfrac{\partial f}{\partial
I}(S_{0},0)}{1+\varphi(\mu_{0}+\mu_{1}+\beta)}&0\\\ &\\\
-\dfrac{\varphi\dfrac{\partial f}{\partial
I}(S_{0},0)}{1+\varphi(\mu_{0}+\nu)}&\dfrac{1}{1+\varphi(\mu_{0}+\nu)}\end{pmatrix}$

Following the method in [4], we represent the matrix $J_{0}$ in the form

$J_{0}=\begin{pmatrix}F+T&0\\\ A&C\end{pmatrix},$

where

$F=\dfrac{\varphi\dfrac{\partial f}{\partial
I}(S_{0},0)}{1+\varphi(\mu_{0}+\mu_{1}+\beta)},\quad
T=\dfrac{1}{1+\varphi(\mu_{0}+\mu_{1}+\beta)},\quad
C=\dfrac{1}{1+\varphi(\mu_{0}+\nu},\quad A=-\dfrac{\varphi\dfrac{\partial
f}{\partial I}(S_{0},0)}{1+\varphi(\mu_{0}+\nu)}.$

It is clear that $F$ and $T$ are non-negative, $F+T$ is irreducible, and
matrices $C$ and $T$ satisfy

$\rho(T)<1,\quad\rho(C)<1.$

Hence, the basic reproduction number of the discrete model (21) can be
computed as

$\mathcal{R}_{0}=\rho\big{(}F(1-T)^{-1}\big{)}=\dfrac{1}{\mu_{0}+\mu_{1}+\beta}\dfrac{\partial
f\bigg{(}\dfrac{\Lambda}{\mu_{0}+\nu},\,\,0\bigg{)}}{\partial I},$

which means that the basic reproduction numbers of the models (21) and (3) are
equal.

Similarly to Proposition 2, it is easy to verify that the NSFD model (21) has
a unique DEE point $E^{*}=(S^{*},I^{*},R^{*})$ if and only if
$\mathcal{R}_{0}>1$.

The following result on the local asymptotic stability of the DFE point of the
model (21) is a direct consequence of Theorem 2.1 in [4].

######  Corollary 1 (Local asymptotic stability of the DFE point).

The DFE point $E^{0}$ of the NSFD model (21) is locally asymptotically stable
if $\mathcal{R}_{0}\leq 1$ and unstable if $\mathcal{R}_{0}>1$.

###### Proof.

From Theorem 2.1 in [4], we conclude that the DFE point $(S_{0},I_{0})$ of
(24) is locally asymptotically stable if $\mathcal{R}_{0}<1$ and unstable if
$\mathcal{R}_{0}>1$. Consequently, the local asymptotic stability of $E^{0}$
of (21) is proved. ∎

######  Theorem 5.

The DEE point $E^{*}$ of the NSFD model (21) is locally asymptotically stable
if $\mathcal{R}_{0}>1$.

###### Proof.

The Jacobian matrix of the system (21) at $E^{*}$ is given by

$J(E^{*})=\begin{pmatrix}\dfrac{1-\varphi\dfrac{\partial f}{\partial
S}(S^{*},I^{*})}{1+\varphi(\mu_{0}+\nu)}&-\dfrac{\varphi\dfrac{\partial
f}{\partial I}(S^{*},I^{*})}{1+\varphi(\mu_{0}+\nu)}&0\\\ &\\\
\dfrac{\varphi\dfrac{\partial f}{\partial
S}(S^{*},I^{*})}{1+\varphi(\mu_{0}+\mu_{1}+\beta)}&\dfrac{1+\varphi\dfrac{\partial
f}{\partial I}(S^{*},I^{*})}{1+\varphi(\mu_{0}+\mu_{1}+\beta)}&0\\\ &\\\
\dfrac{\partial R_{n+1}}{\partial S_{n}}(E^{*})&\dfrac{\partial
R_{n+1}}{\partial I_{n}}(E^{*})&\dfrac{1}{1+\varphi\mu_{0}}\end{pmatrix}.$

Hence, one of three eigenvalue of $J(E^{*})$ is
$\lambda_{3}=\dfrac{1}{1+\varphi\mu_{0}}$ and two remaining eigenvalues are
the ones of the sub-matrix

$J^{*}=\begin{pmatrix}\dfrac{1-\varphi\dfrac{\partial f}{\partial
S}(S^{*},I^{*})}{1+\varphi(\mu_{0}+\nu)}&-\dfrac{\varphi\dfrac{\partial
f}{\partial I}(S^{*},I^{*})}{1+\varphi(\mu_{0}+\nu)}\\\ &\\\
\dfrac{\varphi\dfrac{\partial f}{\partial
S}(S^{*},I^{*})}{1+\varphi(\mu_{0}+\mu_{1}+\beta)}&\dfrac{1+\varphi\dfrac{\partial
f}{\partial I}(S^{*},I^{*})}{1+\varphi(\mu_{0}+\mu_{1}+\beta)}\end{pmatrix}.$

It is clear that $|\lambda_{3}|<1$. Denote by $\lambda_{1}$ and $\lambda_{2}$
the eigenvalues of $J^{*}$. We will show that $|\lambda_{1}|<1$ and
$|\lambda_{2}|<1$. Thanks to Theorem 2.10 in [3], it is sufficient to verify
that

$|Tr(J^{*})|<1+\det(J^{*})<2.$ (25)

Indeed, we first have

$\det(J^{*})=\dfrac{1+\varphi\dfrac{\partial f}{\partial
I}(S^{*},I^{*})-\varphi\dfrac{\partial f}{\partial
S}(S^{*},I^{*})}{\big{[}1+\varphi(\mu_{0}+\mu_{1}+\beta)\big{]}\big{[}1+\varphi(\mu_{0}+\nu)\big{]}}.$

So, it follows from (12) that $\det(J^{*})<1$. On the other hand,

$Tr(J^{*})=\dfrac{\bigg{(}1-\varphi\dfrac{\partial f}{\partial
S}(S^{*},I^{*})\bigg{)}\big{[}1+\varphi(\mu_{0}+\mu_{1}+\beta)\big{]}+\bigg{(}1+\varphi\dfrac{\partial
f}{\partial
I}(S^{*},I^{*})\bigg{)}\big{[}1+\varphi(\mu_{0}+\nu)\big{]}}{\big{[}1+\varphi(\mu_{0}+\mu_{1}+\beta)\big{]}\big{[}1+\varphi(\mu_{0}+\nu)\big{]}}.$

Hence, it is easy to verify that

$1+\det(J^{*})-Tr(J^{*})>0,\quad 1+\det(J^{*})+Tr(J^{*})>0,$

or equivalently,

$|Tr(J^{*})|<1+\det(J^{*}).$

Therefore, (25) is verified. This completes the proof. ∎

######  Remark 1.

Summing up the results constructed in this section, we conclude that the NSFD
model (21) is dynamically consistent to the continuous model (3), with respect
to the positivity, boundedness, basic reproduction number and local asymptotic
stability. Also, as will be suggested by numerical examples in the next
section, the GAS of the model (3) is also preserved by the NSFD model (21).

To end this section, we investigate the convergence and error bounds for the
NSFD scheme (21). We recall that a difference scheme is said to be convergent
of order p if the global error $e_{n}$ with $e_{0}=0$ and
$e_{n}:=y_{n}-y(t_{n})$, satisfies (see [6, Section 3.2])

$e_{n}=\mathcal{O}(\Delta t^{p}),\quad n=1,2,3,\ldots.$

######  Theorem 6.

The NSFD scheme (21) is convergent of order $1$.

###### Proof.

Let $\big{(}S(0),I(0),R(0)\big{)}\in\mathbb{R}_{3}^{+}$ be any initial data
for initial value problem (3). From the boundedness of (3) and (21) we set

$\begin{split}&C_{S}:=\sup_{t\geq 0}S(t),\quad C_{I}:=\sup_{t\geq 0}I(t),\quad
C_{R}:=\sup_{t\geq 0}R(t),\quad D_{S}:=\sup_{n\geq 0}\\{S_{n}\\},\quad
D_{I}:=\sup_{n\geq 0}\\{I_{n}\\},\quad D_{R}:=\sup_{n\geq 0}\\{R_{n}\\},\\\
&\Omega^{CD}:=\Big{\\{}(S,I,R)\big{|}0\leq
S\leq\max\\{C_{S},D_{S}\\},\,\,0\leq I\leq\max\\{C_{I},D_{I}\\},\,\,0\leq
R\leq\max\\{C_{R},D_{R}\\}\Big{\\}}.\end{split}$

For each arbitrary but fixed step size $\Delta t^{*}>0$, we denote

$\Omega^{D}:=\Big{\\{}(S,I,R)\big{|}0\leq S\leq D_{S},\,\,0\leq I\leq
D_{I},\,\,0\leq R\leq D_{R}\Big{\\}}.$

Next, we denote by $F_{1}(S,I,R),F_{2}(S,I,R)$ and $F_{3}(S,I,R)$ the 1st, 2nd
and 3rd equations of the right-side functions of the model (3). By the
boundedness of $S(t),I(t),R(t)$ and the continuity of $g_{1},g_{2},g_{3}$, it
is valid to set

$L_{S}:=\sup_{t\geq 0}|S^{\prime\prime}(t)|,\quad L_{I}:=\sup_{t\geq
0}|I^{\prime\prime}(t)|,\quad L_{R}:=\sup_{t\geq 0}|R(t)|,\quad
L:=M_{S}+M_{I}+M_{R}$ (26)

and

$\begin{split}&m_{j}^{S}:=\max_{(S,I,R)\in\Omega^{CS}}\bigg{|}\dfrac{\partial
F_{j}}{\partial S}(S,I,R)\bigg{|},\quad j=1,2,3,\\\
&m_{j}^{I}:=\max_{(S,I,R)\in\Omega^{CS}}\bigg{|}\dfrac{\partial
F_{j}}{\partial I}(S,I,R)\bigg{|},\quad j=1,2,3,\\\
&m_{j}^{R}:=\max_{(S,I,R)\in\Omega^{CS}}\bigg{|}\dfrac{\partial
F_{j}}{\partial R}(S,I,R)\bigg{|},\quad j=1,2,3,\\\
&m:=\max\Big{\\{}\sum_{j=1}^{3}m_{j}^{S},,\,\sum_{j=1}^{3}m_{j}^{I},\,\,\sum_{j=1}^{3}m_{j}^{R}\Big{\\}}.\end{split}$
(27)

Let us denote by $G_{1}(S_{n},I_{n},R_{n},\Delta
t),G_{2}(S_{n},I_{n},R_{n},\Delta t)$ and $G_{3}(S_{n},I_{n},R_{n},\Delta t)$
the 1st, 2nd and 3rd equations of the system (23), respectively. It is easy to
verify that

$\begin{split}&G_{1}(S_{n},I_{n},R_{n},0)\equiv S_{n},\\\
&G_{2}(S_{n},I_{n},R_{n},0)\equiv I_{n},\\\ &G_{3}(S_{n},I_{n},R_{n},0)\equiv
R_{n},\\\ &\dfrac{\partial G_{j}}{\partial\Delta t}(S_{n},I_{n},R_{n})\equiv
F_{j}(S_{n},I_{n},R_{n}),\quad j=1,2,3.\end{split}$ (28)

Assume that $\varphi^{\prime\prime}(\Delta t)$ exists for all $\Delta t\geq
0$. We set

$\begin{split}&M_{1}:=\max_{(S_{n},I_{n},R_{n})\in\Omega^{D}}\bigg{|}\dfrac{\partial^{2}G_{1}}{\partial\Delta
t^{2}}(S_{n},I_{n},R_{n},\Delta t)\bigg{|},\\\
&M_{2}:=\max_{(S_{n},I_{n},R_{n})\in\Omega^{D}}\bigg{|}\dfrac{\partial^{2}G_{2}}{\partial\Delta
t^{2}}(S_{n},I_{n},R_{n},\Delta t)\bigg{|},\\\
&M_{3}:=\max_{(S_{n},I_{n},R_{n})\in\Omega^{D}}\bigg{|}\dfrac{\partial^{2}G_{3}}{\partial\Delta
t^{2}}(S_{n},I_{n},R_{n},\Delta t)\bigg{|},\\\
&M:=M_{1}+M_{2}+M_{3}.\end{split}$ (29)

Now, using the Taylor’s theorem we obtain

$\begin{split}&S(t_{n+1})=S(t_{n})+\Delta
tF_{1}(S(t_{n}),I(t_{n}),R(t_{n}))+\dfrac{\Delta
t^{2}}{2}S^{\prime\prime}(\xi_{S}),\quad t_{n}<\xi_{S}<t_{n+1},\\\
&I(t_{n+1})=I(t_{n})+\Delta tF_{2}(S(t_{n}),I(t_{n}),R(t_{n}))+\dfrac{\Delta
t^{2}}{2}I^{\prime\prime}(\xi_{I}),\quad t_{n}<\xi_{I}<t_{n+1},\\\
&R(t_{n+1})=R(t_{n})+\Delta tF_{3}(S(t_{n}),I(t_{n}),R(t_{n}))+\dfrac{\Delta
t^{2}}{2}S^{\prime\prime}(\xi_{R}),\quad t_{n}<\xi_{R}<t_{n+1},\\\
\end{split}$ (30)

Next, it follows from the Taylor’s theorem and (28) that

$\begin{split}&S_{n+1}=S_{n}+\Delta tF_{1}(S_{n},I_{n},R_{n})+\dfrac{\Delta
t^{2}}{2}\dfrac{\partial^{2}G_{1}}{\partial\Delta
t^{2}}(S_{n},I_{n},R_{n},\Delta t_{S}),\quad 0<\Delta t_{S}<\Delta t,\\\
&I_{n+1}=I_{n}+\Delta tF_{2}(S_{n},I_{n},R_{n})+\dfrac{\Delta
t^{2}}{2}\dfrac{\partial^{2}G_{2}}{\partial\Delta
t^{2}}(S_{n},I_{n},R_{n},\Delta t_{I}),\quad 0<\Delta t_{I}<\Delta t,\\\
&R_{n+1}=R_{n}+\Delta tF_{3}(S_{n},I_{n},R_{n})+\dfrac{\Delta
t^{2}}{2}\dfrac{\partial^{2}G_{3}}{\partial\Delta
t^{2}}(S_{n},I_{n},R_{n},\Delta t_{R}),\quad 0<\Delta t_{R}<\Delta
t.\end{split}$ (31)

Similarly,

$\begin{split}&\Big{|}F_{j}(S(t_{n}),I(T_{n}),R(t_{n}))-F_{j}(S_{n},I_{n},R_{n})\Big{|}\\\
&=\bigg{|}\dfrac{\partial F_{j}}{\partial
S}(\zeta_{j})(S(t_{n})-S_{n})+\dfrac{\partial F_{j}}{\partial
I}(\zeta_{j})(I(t_{n})-I_{n})+\dfrac{\partial F_{j}}{\partial
R}(\zeta_{j})(R(t_{n})-R_{n})\bigg{|}\\\ &\leq\bigg{|}\dfrac{\partial
F_{j}}{\partial S}(\zeta_{j})\bigg{|}|S(t_{n})-S_{n}|+\bigg{|}\dfrac{\partial
F_{j}}{\partial I}(\zeta_{j})\bigg{|}|I(t_{n})-I_{n}|+\bigg{|}\dfrac{\partial
F_{j}}{\partial R}(\zeta_{j})\bigg{|}|R(t_{n})-R_{n}|\\\ &\leq
m_{j}^{S}|S(t_{n})-S_{n}|+m_{j}^{I}|I(t_{n})-I_{n}|+m_{j}^{R}|R(t_{n})-R_{n}|,\end{split}$
(32)

where $\zeta_{j}$ $(j=1,2,3)$ are points between
$(S(t_{n}),I(t_{n}),R(t_{n}))$ and $(S_{n},I_{n},R_{n})$ and
$m_{j}^{S},m_{j}^{I},m_{j}^{R}$ $(j=1,2,3)$ are given by (27).

Let us denote

$e^{S}_{n}=S(t_{n})-S_{n},\quad e^{I}_{n}=I(t_{n})-I_{n},\quad
e^{R}_{n}=R(t_{n})-R_{n},\quad e_{n}=|e_{n}^{S}|+|e_{n}^{I}|+|e_{n}^{R}|.$

Then, from (30)-(32) we obtain

$\begin{split}&|e_{n+1}^{S}|\leq|e_{n}^{S}|+\Delta
t\big{(}m_{1}^{S}|e_{n}^{S}|+m_{1}^{I}|e_{n}^{I}|+m_{1}^{R}|e_{n}^{R}|\big{)}+\dfrac{\Delta
t^{2}}{2}(L_{1}+M_{1}),\\\ &|e_{n+1}^{I}|\leq|e_{n}^{I}|+\Delta
t\big{(}m_{2}^{S}|e_{n}^{S}|+m_{2}^{I}|e_{n}^{I}|+m_{2}^{R}|e_{n}^{R}|\big{)}+\dfrac{\Delta
t^{2}}{2}(L_{2}+M_{2}),\\\ &|e_{n+1}^{R}|\leq|e_{n}^{R}|+\Delta
t\big{(}m_{3}^{S}|e_{n}^{S}|+m_{3}^{I}|e_{n}^{I}|+m_{3}^{R}|e_{n}^{R}|\big{)}+\dfrac{\Delta
t^{2}}{2}(L_{3}+M_{3}),\end{split}$

which implies that

$e_{n+1}\leq(1+m\Delta t)e_{n}+\dfrac{\Delta t^{2}}{2}(L+M),$ (33)

where $L,m$ and $M$ are defined by (26), (27) and (29), respectively.

For convenience, we set $\tau=L+M$. It follows from (33) that

$\begin{split}e_{n+1}&\leq(1+m\Delta t)e_{n}+\tau\dfrac{\Delta t^{2}}{2}\\\
&\leq(1+m\Delta t)\bigg{[}(1+m\Delta t)e_{n-1}+\tau\dfrac{\Delta
t^{2}}{2}\bigg{]}+\tau\dfrac{\Delta t^{2}}{2}=(1+m\Delta
t)^{2}e_{n-1}+\tau\dfrac{\Delta t^{2}}{2}\Big{[}1+(1+m\Delta t)\Big{]}\\\
&\leq\ldots\leq(1+m\Delta t)^{n}e_{0}+\tau\dfrac{\Delta
t^{2}}{2}\sum_{j=0}^{n}(1+m\Delta t)^{j}=\tau\dfrac{\Delta
t^{2}}{2}\dfrac{(1+m\Delta t)^{n+1}-1}{(1+m\Delta t)-1}=\dfrac{\tau\Delta
t}{2m}\big{[}(1+m\Delta t)^{n+1}-1\big{]}.\end{split}$ (34)

Combining (34) and the well-known inequality $e^{x}\geq x+1$ for $x\geq 0$ we
obtain

$e_{n+1}\leq\dfrac{\tau\Delta
t}{2m}\big{[}e^{m(n+1)t_{n}}-1\big{]}=\dfrac{\tau\Delta
t}{2m}\big{[}e^{mt_{n+1}}-1\big{]}.$ (35)

This is the desired conclusion. The proof is complete. ∎

By the estimate (35) we obtain the following result.

######  Corollary 2.

The following estimate holds for the NSFD scheme (21)

$|S(t_{n})-S_{n}|+|I(t_{n})-I_{n}|+|R(t_{n})-R_{n}|\leq\dfrac{\tau\Delta
t}{2m}\big{[}e^{mt_{n}}-1\big{]}$

for $n=0,1,2,\ldots$.

##  4 Numerical experiments

In this section, we report some numerical examples to support the theoretical
results. For this purpose, we consider the model (1) with the incidence
function

$f(S,I)=\dfrac{\alpha SI}{1+aS+bI+cSI},$

where $\alpha,a,b,c>0$. More precisely, the model under consideration is

$\begin{split}&\dfrac{dS}{dt}=\Lambda-\dfrac{\alpha
SI}{1+aS+bI+cSI}-\big{(}\mu_{0}+\nu\big{)}S,\\\ &\dfrac{dI}{dt}=\dfrac{\alpha
SI}{1+aS+bI+cSI}-\big{(}\mu_{0}+\mu_{1}+\beta\big{)}I,\\\
&\dfrac{dR}{dt}=\beta I+\nu S-\mu_{0}R.\end{split}$ (36)

The basic reproduction number of the model (36) is given by

$\mathcal{R}_{0}=\dfrac{\alpha\Lambda}{(\mu_{0}+\mu_{1}+\beta)(\mu_{0}+\nu+a\Lambda)}.$

######  Example 1 (Dynamics of the NSFD scheme and standard ones).

Consider the model (36) with following parameters. In Table 1, the term ”GAS”
stands for the globally asymptotically stable equilibrium point.

Table 1: The parameters used in Example 1 Parameter | Value | Source | Parameter | Value | Source | GAS
---|---|---|---|---|---|---
$\Lambda$ | $0.8$ | Assumed | $\mu_{0}$ | $0.000232$ | [45] | $E^{0}=(649.35,\,0,\,2798.90)$
$\alpha$ | $0.5$ | Assumed | $\nu$ | $0.001$ | Assumed |
$a$ | $0.8$ | Assumed | $\mu_{1}$ | $0.0000547$ | [54] |
$b$ | $0.9$ | Assumed | $\beta$ | $0.8$ | Assumed |
$c$ | $0.95$ | Assumed | $\mathcal{R}_{0}$ | $0.7795$ | Computed |

We use the standard Euler scheme, the second-order Runge-Kutta (RK2) scheme
and the NSFD scheme (21) with $\varphi(\Delta t)=\Delta t$ to solve the model
with initial data $(S(0),I(0),R(0))=(500,100,2000)$. The numerical solutions
generated by these schemes with various step sizes are depicted in Figures
1-9.

From the obtained results we see that the approximations generated by the
Euler and RK2 schemes are negative and unstable. Hence, the positivity,
boundedness and stability of the HBV model are destroyed. The similar result
can be found in previous works on NSFD schemes for differential equations [15,
16, 17, 18, 19, 20, 31, 32]. Conversely, the approximations generated by the
NSFD scheme preserve the essential qualitative properties of the continuous
model regardless of the chosen step sizes. In particular, the results in
Figures 7-9 indicate that the dynamics of the NSFD scheme does not depend on
the step sizes. This can be explained by the theoretical analyses constructed
in Section 3.

Figure 1: The S-component generated by the Euler and NSFD schemes with $\Delta
t=2.5$ after $80$ iterations. Figure 2: The I-component generated by the Euler
and NSFD schemes with $\Delta t=2.5$ after $80$ iterations. Figure 3: The
R-component generated by the Euler and NSFD schemes with $\Delta t=2.5$ after
$80$ iterations. Figure 4: The S-component generated by the RK2 and NSFD
schemes with $\Delta t=2.5$ after $80$ iterations. Figure 5: The I-component
generated by the RK2 and NSFD schemes with $\Delta t=2.5$ after $80$
iterations. Figure 6: The R-component generated by the RK2 and NSFD schemes
with $\Delta t=2.5$ after $80$ iterations. Figure 7: The S-component generated
by the NSFD scheme with $\Delta t=1.0$, $\Delta t=2.0$ and $\Delta t=5.0$.
Figure 8: The I-component generated by the NSFD scheme with $\Delta t=1.0$,
$\Delta t=2.0$ and $\Delta t=5.0$. Figure 9: The R-component generated by the
NSFD scheme with $\Delta t=1.0$, $\Delta t=2.0$ and $\Delta t=5.0$.

######  Example 2 (Dynamics of the HBV model when $\mathcal{R}_{0}<1$).

The aim of this example is to examine the dynamics of the model (3) when
$\mathcal{R}_{0}<1$. For this purpose, consider the model (3) with following
parameters.

Table 2: The parameters used in Example 2 Parameter | Value | Source | Parameter | Value | Source | GAS
---|---|---|---|---|---|---
$\Lambda$ | $0.232$ | [45] | $\mu_{0}$ | $0.000232$ | [45] | $E^{0}=(126.64,\,0,\,873.36)$
$\alpha$ | $0.0009$ | [45] | $\nu$ | $0.0016$ | Assumed |
$a$ | $0.25$ | Assumed | $\mu_{1}$ | $0.0000547$ | [54] |
$b$ | $0.50$ | Assumed | $\beta$ | $0.25$ | Assumed |
$c$ | $0.75$ | Assumed | $\mathcal{R}_{0}$ | $0.0139$ | Computed |

In this case, the DFE point $E^{0}$ is globally asymptotically stable. We use
the NSFD scheme (3) with $\varphi(\Delta t)=\Delta t$ and $\Delta t=10^{-4}$
to simulate the dynamics of the HBV model. The obtained phase spaces are
depicted in Figure 10. In this figure, each blue curve represents a phase
space corresponding to a specific initial data, the red circle indicates the
location of the DFE point and the green arrows describe the evolution of the
HBV model. It is clear that the GAS of $E^{0}$ as well as the dynamical
properties of the HBV model are confirmed.

Figure 10: The phase spaces of the HBV in Example 2.

######  Example 3 (Dynamics of the HBV model when $\mathcal{R}_{0}>1$).

In this example, we investigate the dynamics of the model (3) when
$\mathcal{R}_{0}>1$. For this reason, consider the model (3) with following
parameters.

Table 3: The parameters used in Example 3 Parameter | Value | Source | Parameter | Value | Source | GAS
---|---|---|---|---|---|---
$\Lambda$ | $0.2$ | Assumed | $\mu_{0}$ | $0.000232$ | [45] | $E^{*}=(19.60,\,6.95,\,1004.80)$
$\alpha$ | $0.005$ | Assumed | $\nu$ | $0.0016$ | Assumed |
$a$ | $0.5$ | Assumed | $\mu_{1}$ | $0.0000547$ | [54] |
$b$ | $0.1$ | Assumed | $\beta$ | $0.025$ | Assumed |
$c$ | $0.2$ | Assumed | $\mathcal{R}_{0}$ | $0.0139$ | Computed |

In this case, the DEE point $E^{*}$ is globally asymptotically stable. We use
the NSFD scheme (3) with $\varphi(\Delta t)=t$ and $t=10^{-4}$ to solve the
HBV model. The obtained phase spaces are depicted in Figure 11. Similarly to
Example 2, the GAS of $E^{0}$ as well as the dynamical properties of the model
are shown.

Figure 11: The phase spaces of the HBV in Example 3.

##  5 Conclusions and open problems

In this work, we have studied dynamics of a generalized hepatitis B epidemic
model and its dynamically consistent NSFD model. The positivity, boundedness,
the basic reproduction number and asymptotic stability properties of the model
have been analyzed rigorously. It was proved, by the Lyapunov stability theory
and the Poincare-Bendixson theorem in combination with the Bendixson-Dulac
criterion, that the DFE point is globally asymptotically stable if the basic
reproduction number $\mathcal{R}_{0}\leq 1$ and the DEE point is globally
asymptotically stable whenever $\mathcal{R}_{0}>1$. Besides, the Mickens’
methodology was applied to formulate a dynamically consistent NSFD model for
the continuous model. By rigorously mathematical analyses we have shown that
the constructed NSFD scheme preserves essential mathematical features of the
continuous model for all finite step sizes. Finally, the numerical experiments
are conducted to illustrate the theoretical findings and to demonstrate
advantages of the NSFD scheme over standard ones. The obtained results
provided important improvements for the ones presented in [36] and [59].

Based on the approach in [36], we can also use the optimal control strategy
for the model (3) to eliminate the spreading of the HBV. This is a very
interesting and important problem with many useful applications in real-world
applications.

Although the constructed NSFD scheme is only convergent of order $1$, its main
advantage is that it can preserve the essential properties of the continuous
model for all finite step sizes. On the other hand, it can be easily combined
with extrapolation techniques or variable step-size strategies to improve the
accuracy. The most important thing is that it can be operate well regardless
of step sizes.

In the near future, optimal control strategies and high-order NSFD schemes for
the model (3) and its extensions will be studied. Additionally, the practice
applications of the HBV model (3) will be also considered.

## References

## References

  * [1] O. Adekanye, T. Washington, Nonstandard finite difference scheme for a Tacoma Narrows Bridge model, Applied Mathematical Modelling 62(2018) 223-236.
  * [2] S. Ahmad, M. Rahman, M. Arfan, On the analysis of semi-analytical solutions of Hepatitis B epidemic model under the Caputo-Fabrizio operator, Chaos, Solitons & Fractals 146 (2021) 110892.
  * [3] L. J. S. Allen, An introduction to mathematical biology, Pearson, 2006.
  * [4] L. J. S. Allen, P. van den Driessche, The basic reproduction number in some discrete-time epidemic models, Journal of Difference Equations and Applications 14(2008) 1127-1147.
  * [5] A. J. Arenas, G. Gonzalez-Parra, B. M. Chen-Charpentier, Construction of nonstandard finite difference schemes for the SI and SIR epidemic models of fractional order, Mathematics and Computers in Simulation 121 (2016) 48-63.
  * [6] U. M. Ascher, L. R. Petzold, Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, Society for Industrial and Applied Mathematics, Philadelphia, 1998.
  * [7] S. Bonhoeffer, R. M. May, G. M. Shaw, M. A. Nowak, Virus dynamics and drug therapy. Proc. Natl Acad. Sci. USA 94(1997) 6971-6976.
  * [8] F. Brauer, The Kermack-McKendrick epidemic model revisited, Mathematical Biosciences 198(2005)119-131.
  * [9] F. Brauer, Mathematical epidemiology: Past, present, and future, Infectious Disease Modelling 2(2017)113-127.
  * [10] J. Calatayud, M. Jornet, An improvement of two nonstandard finite difference schemes for two population mathematical models, Journal of Difference Equations and Applications 27(2021)422-430.
  * [11] L. C. Cardoso, R. F. Camargo, F. L. P. dos Santos, J. P. C. D. Santos, Global stability analysis of a fractional differential system in hepatitis B, Chaos, Solitons & Fractals 143(2021) 110619.
  * [12] L. C. Cardoso, F. L. P. Dos Santos, R. F. Camargo, Analysis of fractional-order models for hepatitis B, Computational and Applied Mathematics 37(2018) 457-4586.
  * [13] J. Cresson, F. Pierret, Non standard finite difference scheme preserving dynamical properties, Journal of Computational and Applied Mathematics 303(2016) 15-30.
  * [14] J. Danane, K. Allali, Z. Hammouch, Mathematical analysis of a fractional differential model of HBV infection with antibody immune response, Chaos, Solitons & Fractals 136(2020) 109787.
  * [15] Quang A Dang, Manh Tuan Hoang, Positivity and global stability preserving NSFD schemes for a mixing propagation model of computer viruses, Journal of Computational and Applied Mathematics 374(2020) 112753.
  * [16] Quang A Dang, Manh Tuan Hoang, Nonstandard finite difference schemes for a general predator-prey system, Journal of Computational Science 36(2019) 101015.
  * [17] Quang A Dang, Manh Tuan Hoang, Dynamically consistent discrete metapopulation model, Journal of Difference Equations and Applications 22 (2016) 1325-1349.
  * [18] Quang A Dang, Manh Tuan Hoang, Lyapunov direct method for investigating stability of nonstandard finite difference schemes for metapopulation models, Journal of Difference Equations and Applications 24(2019) 15-47.
  * [19] Quang A Dang, Manh Tuan Hoang, Complete global stability of a metapopulation model and its dynamically consistent discrete models, Qualitative theory of dynamical systems 18 (2019) 461-475.
  * [20] Quang A Dang, Manh Tuan Hoang, Positive and elementary stable explicit nonstandard Runge-Kutta methods for a class of autonomous dynamical systems, International Journal of Computer Mathematics 97 (2020) 2036-2054.
  * [21] A. Din, Y. Li, Q. Liu, Viral dynamics and control of hepatitis B virus (HBV) using an epidemic model, Alexandria Engineering Journal 59(2020) 667-679.
  * [22] A. Din, Y. Li, A. Yusuf, Delayed hepatitis B epidemic model with stochastic analysis, Chaos, Solitons & Fractals 146 (2021) 110839.
  * [23] W. D. Qin, Q. Ma, Z. Y. Man, X. H. Ding, A boundedness and monotonicity preserving method for a generalized population model, Journal of Difference Equations and Applications 26(2020)1347-1368.
  * [24] F. Gao, X. Li, W. Li, X. Zhou, Stability analysis of a fractional-order novel hepatitis B virus model with immune delay based on Caputo-Fabrizio derivative, Chaos, Solitons & Fractals 142(2021) 110436.
  * [25] S. M. Garba, A. B. Gumel, A. S. Hassan, J. M. -S. Lubuma, Switching from exact scheme to nonstandard finite difference scheme for linear delay differential equation, Applied Mathematics and Computation 258(2015) 388-403.
  * [26] M. Gupta, J. M. Slezak, F. Alalhareth, S. Roy, H. V. Kojouharov, Second-order Nonstandard Explicit Euler Method, AIP Conference Proceedings 2302(2020) 110003.
  * [27] N. C. Grassly, C. Fraser, Mathematical models of infectious disease transmission, Nature Reviews Microbiology 6(2008) 477-487.
  * [28] H. W. Hethcote, The mathematics of infectious diseases. SIAM Rev. 42((2000) 599-653.
  * [29] M. T. Hoang, O. F. Egbelowo O.F, Dynamics of a Fractional-Order Hepatitis B Epidemic Model and Its Solutions by Nonstandard Numerical Schemes. In: Hattaf K., Dutta H. (eds) Mathematical Modelling and Analysis of Infectious Diseases. Studies in Systems, Decision and Control, vol 302. Springer, Cham. https://doi.org/10.1007/978-3-030-49896-2_5.
  * [30] M. T. Hoang, O.F. Egbelowo, On the global asymptotic stability of a hepatitis B epidemic model and its solutions by nonstandard numerical schemes, Boletín de la Sociedad Matemática Mexicana 26(2020) 1113-1134.
  * [31] Manh Tuan Hoang, Reliable approximations for a hepatitis B virus model by nonstandard numerical schemes, Mathematics and Computers in Simulation 193(2022) 32-56.
  * [32] Manh Tuan Hoang, A. M. Nagy, Uniform asymptotic stability of a Logistic model with feedback control of fractional order and nonstandard finite difference schemes, Chaos, Solitons & Fractals 123(2019) 24-34.
  * [33] Z. Horváth, Positivity of Runge-Kutta and diagonally split Runge-Kutta methods, Applied Numerical Mathematics 28(1998) 309-326.
  * [34] P. T. Karaji, N. Nyamoradi, Analysis of a fractional SIR model with General incidence function, Applied Mathematics Letters 108(2020) 106499.
  * [35] W. O. Kermack, A. G. McKendrick, A Contribution to the Mathematical Theory of Epidemics, Proceedings of the Royal Society of London, Series A 115(1927) 700-721.
  * [36] T. Khan, Z. Ullah, N. Ali, G. Zaman, Modeling and control of the hepatitis B virus spreading using an epidemic model, Chaos, Solitons and Fractals 124 (2019) 1-9.
  * [37] A. Khan, G. Hussain, M. Inc, G. Zaman, Existence, uniqueness, and stability of fractional hepatitis B epidemic model, Chaos 30(2020) 103104.
  * [38] T. Khan, A. Khan, G. Zaman, The extinction and persistence of the stochastic hepatitis B epidemic model, Chaos, Solitons & Fractals 108(2018) 123-128.
  * [39] H. V. Kojouharov, S. Roy, M. Gupta, F. Alalhareth, J. M. Slezak, A second-order modified nonstandard theta method for one-dimensional autonomous differential equations, Applied Mathematics Letters 112(2021) 106775.
  * [40] J. La Salle, S. Lefschetz, Stability by Liapunov’s Direct Method, Academic Press, New York, 1961.
  * [41] M. Y. Li, An Introduction to Mathematical Modeling of Infectious Diseases, Springer International Publishing AG, 2018.
  * [42] A. M. Lyapunov, The general problem of the stability of motion, International Journal of Control, Taylor & Francis, 1992.
  * [43] K. Manna, S. P. Chakrabarty, Global stability of one and two discrete delay models for chronic hepatitis B infection with HBV DNA-containing capsids, Computational and Applied Mathematics 36(2017) 525-536.
  * [44] M. Martcheva, An Introduction to Mathematical Epidemiology, Springer, New York, 2015.
  * [45] N. K. Martin, P. Vickerman, M. Hickman, Mathematical modelling of hepatitis C treatment for injecting drug users, Journal of Theoretical Biology 274(2011) 58-66.
  * [46] C. C. McCluskeya, Y. Yang, Global stability of a diffusive virus dynamics model with general incidence function and time delay, Nonlinear Analysis: Real World Applications 25(2015) 64-78.
  * [47] A. McNabb, Comparison theorems for differential equations, Journal of Mathematical Analysis and Applications 119(1986) 417-428.
  * [48] R. E. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, 1993.
  * [49] R. E. Mickens, Applications of Nonstandard Finite Difference Schemes, World Scientific, 2000.
  * [50] R. E. Mickens, Advances in the Applications of Nonstandard Finite Difference Schemes, World Scientific, 2005.
  * [51] R. E. Mickens, Nonstandard Finite Difference Schemes for Differential Equations, Journal of Difference Equations and Applications 8(2002) 823-847.
  * [52] R. E. Mickens, Nonstandard Finite Difference Schemes: Methodology and Applications, World Scientific, 2020.
  * [53] R. E. Mickens, I. H. Herron, Approximate rational solutions to the Thomas-Fermi equation based on dynamic consistency, Applied Mathematics Letters 116(2021)106994.
  * [54] MMWR, Progress in hepatitis B prevention through universal infant vaccination China, 1997-2006, Morbidity and mortality weekly report 2007(18) 441-5 .
  * [55] M. A. Nowak, R. M. May, Virus Dynamics: Mathematical Principles of Immunology and Virology, New York: Oxford University Press 2000.
  * [56] K. C. Patidar, On the use of nonstandard finite difference methods, Journal of Difference Equations and Applications 11(2005) 735-758.
  * [57] K. C. Patidar, Nonstandard finite difference methods: recent trends and further developments, Journal of Difference Equations and Applications 22(2016) 817-849.
  * [58] S. A. A. Shah, M. A. Khan, M. Farooq, S. Ullah, E. O.Alzahrani, A fractional order model for Hepatitis B virus with treatment via Atangana-Baleanu derivative, Physica A: Statistical Mechanics and its Applications 538 (2020) 122636.
  * [59] A. Suryanto, I. Darti, On the nonstandard numerical discretization of SIR epidemic model with a saturated incidence rate and vaccination, AIMS Mathematics 6(2020) 141-155.
  * [60] Y. Tian, X. Liu, Global dynamics of a virus dynamical model with general incidence rate and cure rate, Nonlinear Analysis: Real World Applications 16(2014) 17-26.
  * [61] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathematical Biosciences 180(2002) 29-48.
  * [62] D. T. Wood, H. V. Kojouharov, A class of nonstandard numerical methods for autonomous dynamical systems, Applied Mathematics Letters 50(2015) 78-82.
  * [63] D. T. Wood, H. V. Kojouharov, D. T. Dimitrov, Universal approaches to approximate biological systems with nonstandard finite difference methods, Mathematics and Computers in Simulation 133(2017) 337-350.