Nonstandard Finite Difference Schemes for the Study of the Dynamics of the Babesiosis Disease

Abstract

In this paper, a discrete-time model for Babesiosis disease, given by means of nonstandard finite difference (NSFD) schemes, is first provided and analyzed. Mathematical analyses show that the provided NSFD schemes preserve the essential (qualitative) dynamical properties of the continuous-time model, namely, positivity and boundedness of the solutions, equilibria, and their stability properties. In particular, the global stability of the disease free equilibrium point is proved by using an appropriate Lyapunov function. As a relevant consequence, we get the dynamic consistency of NSFD schemes in relation to the continuous-time model. Numerical simulations are presented to support the validity of the established theoretical results.

1. Introduction

The bovine babesiosis, caused by Babesia bovis and Babesia bigemina, is one of the most important vector-transmitted diseases. It is transmitted by the sting of ticks as the principal vector [1,2]. Babesiosis has an important economic impact in the livestock sector of tropical regions. Mathematical models for Babesiosis disease play an important role in both theory and practice. Therefore, these models have attracted the attention of many mathematicians, biologists, ecologists, and epidemiologists on several aspects, including the qualitative study (see in [3,4,5,6,7,8,9] and references therein) and the construction of adequate discretizations [10].

In this paper, we consider the well-known continuous-time mathematical model for Babesiosis disease proposed in [3], in order to construct nonstandard finite difference (NSFD) schemes, preserving the essential properties (i.e., positivity, boundedness, and stability) of the continuous-time model. It should be emphasized that the transformation of continuous models into discrete ones with the preservation of the essential properties is very important but not simple work. The most successful approach to the problem is the use of NSFD schemes. So far, NSFD schemes have been one of the most effective tools for this problem, as the majority of standard finite difference schemes can change and not preserve the properties of the continuous model [11] for any grid size. Notable results on NSFD for systems of ODEs can be found in [11,12,13,14,15,16,17,18,19,20]. There, the considered systems are the models of essential phenomena and processes arising in applied fields. Recently, some results on NSFD schemes for some ordinary and fractional differential equations have been obtained [21,22,23,24,25,26,27].

Returning to the work [3], notice that some stability properties of the model were established either theoretically or numerically. In [10], a discrete-time model for Babesiosis, which is in essence the Euler scheme with the unit step size, was studied and similar results of stability as in the continuous-time model were obtained either by rigorous proof or numerical simulations. Nevertheless, in the discrete case, the additional assumption of some parametric constrains was needed. Motivated by this work, in this paper, we construct discrete models for system (2) with all the stability properties established theoretically and illustrated by numerical examples. The numerical simulations also confirm the advantages of NSFD schemes over standard finite difference (SFD) schemes for large grid sizes.

By this occasion, we remark that a big challenge in the construction of NSFD schemes for continuous-time epidemic models is that the disease-free equilibrium point of the models is not only locally stable but also globally stable. Furthermore, when the reproduction number ${\mathcal{R}}_0=1$, it becomes a non-hyperbolic equilibrium point. Many continuous-time models considered before (see, e.g., in [11,13,14,15,16,20]) only deal with the local stability of hyperbolic equilibrium points. Here, by using an appropriate Lyapunov function, we prove that the disease-free equilibrium point is globally asymptotically stable if ${\mathcal{R}}_0\le 1$. As the main conclusion of this study, we get that NSFD schemes are dynamically consistent with the continuous model. On the other hand, our results in this paper constitute a generalization of those obtained in the recent work [10].

This paper is organized as follows. The mathematical model and its properties are recalled briefly in Section 2. Section 3 is devoted to the construction of NSFD schemes. Section 4 discusses the numerical simulations. Finally, there are some conclusions and discussions.

2. Mathematical Model

First, we consider the mathematical model for Babesiosis disease constructed in [3]. The dynamic transmission of Babesiosis disease for bovine and tick populations can be modeled by the following system of nonlinear first order differential equations,

$$\left\{\begin{array}{cc}{\overline{S^{\prime }}}_B\left(t\right)\hfill & =\left({\mu }_B+{\alpha }_B\right){\overline{C}}_B\left(t\right)-{\beta }_B{\overline{S}}_B\left(t\right){\displaystyle \frac{{\overline{I}}_T\left(t\right)}{N_T\left(t\right)}},\hfill \\ {\overline{I^{\prime }}}_B\left(t\right)\hfill & ={\beta }_B{\overline{S}}_B\left(t\right){\displaystyle \frac{{\overline{I}}_T\left(t\right)}{N_T\left(t\right)}}-{\lambda }_B{\overline{I}}_B\left(t\right),\hfill \\ {\overline{C^{\prime }}}_B\left(t\right)\hfill & ={\lambda }_B{\overline{I}}_B\left(t\right)-\left({\mu }_B+{\alpha }_B\right){\overline{C}}_B\left(t\right),\hfill \\ {\overline{S^{\prime }}}_T\left(t\right)\hfill & ={\mu }_{T}p{\overline{I}}_T-{\overline{S}}_T\left(t\right){\displaystyle \frac{{\overline{I}}_B\left(t\right)}{N_B\left(t\right)}},\hfill \\ {\overline{I^{\prime }}}_T\left(t\right)\hfill & ={\overline{S}}_T\left(t\right){\displaystyle \frac{{\overline{I}}_B\left(t\right)}{N_B\left(t\right)}}-{\mu }_{T}p{\overline{I}}_T,\hfill \end{array}\right.$$

(1)

where the total population of bovine $N_B\left(t\right)$ is divided into three subpopulations: bovines which may become infected (Susceptible ${\overline{S}}_B\left(t\right)$), bovines infected by the Babesia parasite (Infected ${\overline{I}}_B\left(t\right)$), and bovines which have been treated for the Babesiosis (Controlled ${\overline{C}}_B\left(t\right)$). That is, as in [3], $N_B=\overline{S_B}+\overline{I_B}+\overline{C_B}$. The total population of ticks $N_T\left(t\right)$ is divided into two subpopulations: ticks which may become infected ${\overline{S}}_T\left(t\right)$ and ticks infected by the Babesia parasite ${\overline{I}}_T\left(t\right)$. As in [3], the parameter ${\mu }_B$ represents the birth rate of the bovines and it is assumed equal to their natural death; ${\mu }_T$ is the birth rate of the ticks and it is also assumed equal to their death rate; any susceptible bovine can become infected due to a sting of an infected tick at a rate ${\beta }_B$; any susceptible tick can become infected when it stings an infected bovine, at a rate ${\beta }_T$; p is the probability for a susceptible tick to be born from an infected one; ${\lambda }_B$ represents the fraction of the infected bovines which are controlled, that is, treated against Babesia, while ${\alpha }_B$ represents the fraction of the controlled bovines which may return to susceptible state. All these parameters in the model are positive and not exceeding 1, as this is biologically logical. More details on this model can be found in [3].

Below, we recall the main results of that work. In particular, after normalization of the variables, it was shown that system (1) is equivalent to subsystem

$$\left\{\begin{array}{cc}{S}_B^{\prime }\left(t\right)\hfill & =\left({\mu }_B+{\alpha }_B\right)\left(1-S_B\left(t\right)-I_B\left(t\right)\right)-{\beta }_B{S}_B\left(t\right)I_T\left(t\right),\hfill \\ I_B^{\prime }\left(t\right)\hfill & ={\beta }_B{S}_B\left(t\right)I_T\left(t\right)-{\lambda }_B{I}_B\left(t\right),\hfill \\ I_T^{\prime }\left(t\right)\hfill & ={\beta }_T\left(1-I_T\left(t\right)\right)I_B\left(t\right)-{\mu }_{T}p{I}_T\left(t\right)\hfill \end{array}\right.$$

(2)

and the region

$$\mathsf{\Omega }=\left\{\left(S_B,I_B,I_T\right)\in {\mathbb{R}}_{+}^3|0\le S_B+I_B\le 1,0\le I_T\le 1\right\},$$

(3)

is a positive invariant set for the system (2). For this system, the threshold parameter is given by

$${\mathcal{R}}_0={\displaystyle \frac{{\beta }_B{\beta }_T}{{\lambda }_B{\mu }_{T}p}}.$$

(4)

Concerning system (2), the following results are obtained.

(i)

System (2) has the disease-free equilibrium point $F_1^{*}=(S_{B1}^{*},I_{B1}^{*},I_{T1}^{*})=(1,0,0)$ for all values of the parameters in this system, whereas, only if ${\mathcal{R}}_0>1$, there is a (unique) endemic equilibrium point $F_2^{*}=(S_{B2}^{*},I_{B2}^{*},I_{T2}^{*})$ in the interior of $\mathsf{\Omega }$ given by

$$\left\{\begin{array}{cc}{S}_{B2}^{*}\hfill & ={\displaystyle \frac{{\beta }_T{\lambda }_B({\alpha }_B+{\mu }_B)+p{\lambda }_B({\alpha }_B+{\lambda }_B+{\mu }_B){\mu }_T}{{\beta }_T\left[{\alpha }_B({\beta }_B+{\lambda }_B)+{\lambda }_B+{\mu }_B+{\beta }_B({\lambda }_B+{\mu }_B)\right]}},\hfill \\ I_{B2}^{*}\hfill & ={\displaystyle \frac{({\alpha }_B+{\mu }_B)({\beta }_B{\beta }_T-p{\lambda }_B{\mu }_T)}{{\beta }_T\left[{\alpha }_B({\beta }_B+{\lambda }_B)+{\lambda }_B+{\mu }_B+{\beta }_B({\lambda }_B+{\mu }_B)\right]}},\hfill \\ I_{T2}^{*}\hfill & ={\displaystyle \frac{({\alpha }_B+{\mu }_B)({\beta }_B{\beta }_T-p{\lambda }_B{\mu }_T)}{{\beta }_T{\beta }_B({\alpha }_B+{\mu }_B)+p{\beta }_B({\alpha }_B+{\lambda }_B+{\mu }_B){\mu }_T}}.\hfill \end{array}\right.$$

(5)

(ii)

If ${\mathcal{R}}_0\le 1$, then the disease-free point $F_1^{*}$ is globally asymptotically stable; otherwise, the disease-free point $F_1^{*}$ is unstable.

(iii)

If ${\mathcal{R}}_0>1$, then the endemic point $F_2^{*}$ is shown to be locally asymptotically stable by numerical simulations.

In [10], the authors proposed a discrete-time version of the continuous model (2) and obtained similar results of stability as for the continuous case. However, the additional assumption of some parametric constraints was required.

3. Nonstandard Finite Difference Schemes for System (2)

Our main goal is to construct NSFD schemes preserving the essential properties of the model (2). Let N be a positive integer and $[0,T]$ be a finite interval. Let us denote by $h=\mathsf{\Delta }t=T/N$ the time step size of the discretization:

$$0=t_0<t_1<\dots <t_N=T=Nh,$$

and let $t_k=kh$ for $k=0,1,\dots ,N$. We now denote by $S_B^k,I_B^k$ and $I_T^k$ the approximated values of $S_B\left(t_k\right),I_B\left(t_k\right)$ and $I_T\left(t_k\right)$, respectively. Following the Mickens’ methodology, we propose NSFD schemes of the form

$$\left\{\begin{array}{cc}{\displaystyle \frac{S_B^{k+1}-S_B^k}{\phi \left(h\right)}}\hfill & =\left({\mu }_B+{\alpha }_B\right)\left(1-S_B^k-I_B^k\right)-{\beta }_B{S}_B^k{I}_T^k,\hfill \\ {\displaystyle \frac{I_B^{k+1}-I_B^k}{\phi \left(h\right)}}\hfill & ={\beta }_B{S}_B^k{I}_T^k-{\lambda }_B{I}_B^k,\hfill \\ {\displaystyle \frac{I_T^{k+1}-I_T^k}{\phi \left(h\right)}}\hfill & ={\beta }_T\left(1-I_T^k\right)I_B^k-{\mu }_{T}p{I}_T^k,\hfill \end{array}\right.$$

(6)

where $\phi \left(h\right)=h+\mathcal{O}\left(h^2\right)$ as $h\to 0$. The conditions on $\phi \left(h\right)$ will be determined so that the properties of system (2) are preserved. For convenience, h is omitted in some presentations.

Theorem 1.

The set Ω defined by (3) is a positive invariant set of discrete model (6) if

$$\phi \left(h\right)<min\left\{{\displaystyle \frac{1}{{\beta }_B}},{\displaystyle \frac{1}{{\lambda }_B}},{\displaystyle \frac{1}{{\mu }_{T}p}},\frac{1}{{\beta }_T},{\displaystyle \frac{1}{{\mu }_B+{\alpha }_B}}\right\},\forall h>0.$$

(7)

Proof. 

The theorem is proved by induction. First, rewrite system (6) in the form

$$\left\{\begin{array}{cc}{S}_B^{k+1}\hfill & =S_B^k+\phi \left({\mu }_B+{\alpha }_B\right)\left(1-S_B^k-I_B^k\right)-\phi {\beta }_B{I}_T^k{S}_B^k,\hfill \\ I_B^{k+1}\hfill & =I_B^k+\phi {\beta }_B{S}_B^k{I}_T^k-\phi {\lambda }_B{I}_B^k,\hfill \\ I_T^{k+1}\hfill & =I_T^k+\phi {\beta }_T\left(1-I_T^k\right)I_B^k-\phi {\mu }_{T}p{I}_T^k.\hfill \end{array}\right.$$

(8)

Clearly, if $(S_B^k,I_B^k,I_T^k)\in \mathsf{\Omega }$ and condition (7) holds, then $S_B^{k+1}\ge 0$, $I_B^{k+1}\ge 0$ and $I_T^{k+1}\ge 0$. On the other hand, given $(S_B^k,I_B^k,I_T^k)\in \mathsf{\Omega }$, from the third equation of (8), we have

$$I_T^{k+1}\le I_T^k+\phi {\beta }_T\left(1-I_T^k\right)=(1-\phi {\beta }_T)I_T^k+\phi {\beta }_T\le 1-\phi {\beta }_T+\phi {\beta }_T=1.$$

Adding the first and second equations of (8), we obtain

$$\begin{array}{cc}\hfill S_B^{k+1}+I_B^{k+1}& =S_B^k+I_B^k+\phi \left({\mu }_B+{\alpha }_B\right)\left(1-S_B^k-I_B^k\right)-\phi {\lambda }_B{I}_B^k\le S_B^k+I_B^k+\phi \left({\mu }_B+{\alpha }_B\right)\left(1-S_B^k-I_B^k\right)\hfill \\ \hfill & =\left[1-\phi ({\mu }_B+{\alpha }_B)\right](S_B^k+I_B^k)+\phi ({\mu }_B+{\alpha }_B)\le 1-\phi ({\mu }_B+{\alpha }_B)+\phi ({\mu }_B+{\alpha }_B)=1.\hfill \end{array}$$

Thus, the proof is complete. □

Similar to Proposition 1 in [3], we obtain the following result.

Corollary 1.

Model (6) has the disease-free equilibrium point$F_1^{*}=(S_B^{*},I_B^{*},I_T^{*})=(1,0,0)$for all the values of the parameters in this model, whereas, only if${\mathcal{R}}_0>1$, there is a (unique) endemic equilibrium point$F_2^{*}=(S_B^{*},I_B^{*},I_T^{*})$in the interior of Ω, where the endemic point is given by (5).

Theorem 2.

Consider the model (6) under the assumption of Theorem 1.

(i) 

If${\mathcal{R}}_0\le 1$, then the disease-free equilibrium point$F_1^{*}$is globally asymptotically stable.

(ii) 

If${\mathcal{R}}_0>1$, then the disease-free equilibrium point$F_1^{*}$is unstable.

Proof. 

We shall distinguish two parts as in the statement of the theorems.

(i) 

We will use an extension for the discrete case (see in [28], Theorem 3.3) of the Lyapunov stability theorem [29] to prove this part. For this purpose, consider a function $V:\mathsf{\Omega }\to {\mathbb{R}}_{+}$ defined by

$$V(S_B^k,I_B^k,I_T^k):={\beta }_T{I}_B^k+{\lambda }_B{I}_T^k.$$

Clearly, V is continuous, $V\left(S_B^k,I_B^k,I_T^k\right)\ge 0$ for all $\left(S_B^k,I_B^k,I_T^k\right)\in \mathsf{\Omega }$, and $V\left(F_1^{*}\right)=0$.

From (6), we have

$$\begin{array}{cc}\hfill \mathsf{\Delta }V(S_B^k,I_B^k,I_T^k)& :=V(S_B^{k+1},I_B^{k+1},I_T^{k+1})-V(S_B^k,I_B^k,I_T^k)={\beta }_T(I_B^{k+1}-I_B^k)+{\lambda }_B(I_T^{k+1}-I_T^k),\hfill \\ \hfill & =\phi {\beta }_T\left({\beta }_B{S}_B^k{I}_T^k-{\lambda }_B{I}_B^k\right)+\phi {\lambda }_B\left[{\beta }_T\left(1-I_T^k\right)I_B^k-{\mu }_{T}p{I}_T^k\right]\hfill \\ \hfill & \le \phi {\beta }_T\left({\beta }_B{I}_T^k-{\lambda }_B{I}_B^k\right)+\phi {\lambda }_B\left[{\beta }_T\left(1-I_T^k\right)I_B^k-{\mu }_{T}p{I}_T^k\right]\hfill \\ \hfill & =\phi \left[({\beta }_T{\beta }_B-{\lambda }_B{\mu }_{T}p)I_T^k-{\lambda }_B{\beta }_T{I}_B^k{I}_T^k\right]\hfill \\ \hfill & \le \phi ({\beta }_T{\beta }_B-{\lambda }_B{\mu }_{T}p)I_T^k=\phi {\lambda }_B{\mu }_{T}p({\mathcal{R}}_0-1)I_T^k,\hfill \end{array}$$

(9)

which implies that $\mathsf{\Delta }V\le 0$ for all $\left(S_B^k,I_B^k,I_T^k\right)\in \mathsf{\Omega }$.

Let $G^{*}$ be the largest positively invariant set contained in

$$G:=\left\{\left(S_B^k,I_B^k,I_T^k\right)\in \mathsf{\Omega }|\mathsf{\Delta }V=0\right\}.$$

Then, by using (9) we have that

$$G^{*}=\left\{\begin{array}{c}\left\{F_1^{*}\right\}\mathrm{if}{\mathcal{R}}_0<1,\hfill \\ \\ \left\{(0,0,I_T^k)|I_T^k\ge 0\right\}\mathrm{if}{\mathcal{R}}_0=1.\hfill \end{array}\right.$$

Consequently, it is easy to verify that $F_1^{*}$ is $G^{*}$-globally asymptotically stable if ${\mathcal{R}}_0\le 1$.

As all solutions of (6) are bounded, by (see in [28], Theorem 3.3) we deduce that $F_1^{*}$ is globally asymptotically stable if ${\mathcal{R}}_0\le 1$.

(ii) 

Computing the Jacobian matrix of system (6) evaluated at the disease free point, one obtains

$$J(1,0,0)=\left(\begin{array}{ccc}1-\phi ({\alpha }_B+{\mu }_B)& -\phi ({\alpha }_B+{\mu }_B)& -\phi {\beta }_B,\\ & & \\ 0& 1-\phi {\lambda }_B& \phi {\beta }_B\\ & & \\ 0& \phi {\beta }_T& 1-\phi {\mu }_{T}p\end{array}\right).$$

Consequently, the eigenvalues of $J\left(F_1^{*}\right)$ are ${\mathsf{\Lambda }}_1=1-\phi ({\alpha }_B+{\mu }_B)$ and ${\mathsf{\Lambda }}_2,{\mathsf{\Lambda }}_3$, where ${\mathsf{\Lambda }}_2,{\mathsf{\Lambda }}_3$ are the eigenvalues of

$$J_1=\left(\begin{array}{cc}1-\phi {\lambda }_B& \phi {\beta }_B\\ & & \\ \phi {\beta }_T& 1-\phi {\mu }_{T}p\end{array}\right).$$

We have

$$det\left(J_1\right)=1-\phi ({\lambda }_B+{\mu }_{T}p)-{\phi }^2{\lambda }_B{\mu }_{T}p({\mathcal{R}}_0-1),Tr\left(J_1\right)=2-\phi ({\lambda }_B+{\mu }_{T}p).$$

Thus, if ${\mathcal{R}}_0>1$, then $1-Tr\left(J_1\right)+det\left(J_1\right)=-{\phi }^2{\lambda }_B{\mu }_{T}p({\mathcal{R}}_0-1)<0$. By Theorem 1.3.7 in [30] and Theorem 2.10 in [31], we can conclude that if ${\mathcal{R}}_0>1$ then the disease-free equilibrium point $F_1^{*}$ is unstable.

Thus, the theorem is proved. □

Next, we shall investigate the local stability of system (6) in the case ${\mathcal{R}}_0>1$. Denote by $\sigma \left(J\right)$ the set of eigenvalues of the Jacobian matrix J of system (2) at $F_2^{*}$. Based on mathematical analyses in [3], we have if ${\mathcal{R}}_0>1$ then $F_2^{*}$ is locally asymptotically stable and the eigenvalues of J lie within the left half of the complex plane, i.e.,

$$Re\left({\mathsf{\Lambda }}_i\right)<0,{\mathsf{\Lambda }}_i\in \sigma \left(J\right),i=1,2,3.$$

Now we denote by K the Jacobian matrix of system (6) at $F_2^{*}$. Then $K=I+\phi J$ (see [13,14]), where I is $3\times 3$ unit matrix. Therefore, the eigenvalues ${\mathsf{\Lambda }}_i$ of J correspond to the eigenvalues ${\mathsf{\Psi }}_i$ of K determined by

$${\mathsf{\Psi }}_i=1+\phi {\mathsf{\Lambda }}_i,i=1,2,3.$$

By Lyapunov indirect method for discrete dynamical systems (see ([30], Theorem 1.3.7) and [31]), $F_2^{*}$ is locally asymptotically stable of system (6) if and only if

$$|{\mathsf{\Psi }}_i|<1,i=1,2,3.$$

This is equivalent to

$$|{\mathsf{\Psi }}_i{|}^2={|1+\phi {\mathsf{\Lambda }}_i|}^2={\phi }^{2}I{m}^2\left({\mathsf{\Lambda }}_i\right)+1+2\phi Re\left({\mathsf{\Lambda }}_i\right)+{\phi }^{2}R{e}^2\left({\mathsf{\Lambda }}_i\right)<1,i=1,2,3.$$

The above inequality is satisfied if

$$\phi <-{\displaystyle \frac{2Re\left({\mathsf{\Lambda }}_i\right)}{|{\mathsf{\Lambda }}_i{|}^2}},i=1,2,3.$$

Conversely, if ${\mathcal{R}}_0\le 1$, then $F_2^{*}$ has no biological sense. Thus, it has no biological sense to study its stability.

From the above results, we have the following theorem.

Theorem 3.

Consider system (6) when${\mathcal{R}}_0>1$. Suppose${\mathsf{\Lambda }}_1,{\mathsf{\Lambda }}_2,{\mathsf{\Lambda }}_3$are the eigenvalues of the Jacobian matrix of system (2) at$F_2^{*}$. Set

$${\tau }^{*}:=\underset{i=1,2,3}{min}\left\{-{\displaystyle \frac{2Re\left({\mathsf{\Lambda }}_i\right)}{|{\mathsf{\Lambda }}_i{|}^2}}\right\}.$$

(10)

Then,$F_2^{*}$is locally asymptotically stable if

$$0<\phi \left(h\right)<{\tau }^{*},\forall h>0.$$

(11)

Summing up the results of this section, we obtain the following theorem.

Theorem 4 (Dynamically consistent discrete models).

Consider NSFD scheme (6).

(i) 

In the case${\mathcal{R}}_0\le 1$, scheme (6) preserves positivity, boundedness, and global stability of$F_1^{*}$of system (2) if

$$0<\phi \left(h\right)<{\phi }^{*}:=min\left\{{\displaystyle \frac{1}{{\beta }_B}},{\displaystyle \frac{1}{{\lambda }_B}},{\displaystyle \frac{1}{{\mu }_{T}p}},\frac{1}{{\beta }_T},{\displaystyle \frac{1}{{\mu }_B+{\alpha }_B}}\right\},\forall h>0.$$

(ii) 

In the case${\mathcal{R}}_0>1$, scheme (6) preserves positivity, boundedness and local stability of$F_2^{*}$and unstability of$F_1^{*}$if

$$0<\phi \left(h\right)<{\phi }^{*}:=min\left\{{\displaystyle \frac{1}{{\beta }_B}},{\displaystyle \frac{1}{{\lambda }_B}},{\displaystyle \frac{1}{{\mu }_{T}p}},\frac{1}{{\beta }_T},{\displaystyle \frac{1}{{\mu }_B+{\alpha }_B}},{\tau }^{*}\right\},\forall h>0,$$

where${\tau }^{*}$is given by (10).

Remark 1.

Concerning Theorem 2, we have to underline three important observations:

• Part (i) of Theorem 2 is only appropriate when $\mathsf{\Omega }$ is a positively invariant set of (6).
• There are many ways for selecting the function $\phi \left(h\right)$ satisfying the conditions of Theorem 4, for example (see in [11,13,14,15,16]),

  $$\phi \left(h\right)={\displaystyle \frac{1-e^{-T^{*}h}}{T^{*}}},T^{*}>{\displaystyle \frac{1}{{\phi }^{*}}}.$$
• In the case of $\phi \left(h\right)=1$, then scheme (6) becomes the discrete model which was proposed and studied in [10]. This proves that NSFD schemes (6) constitute a generalization of the discrete-time model in [10].

4. Numerical Simulations

We present some numerical simulations to illustrate the obtained theoretical results. From the numerical simulations, it will be shown that standard finite difference schemes (SFDS) may not preserve essential properties of the continuous-time model for any finite step size. Meanwhile, NSFD schemes preserve essential properties of continuous model for any finite step size.

Example 1 (Dynamics of standard finite difference schemes).

Consider system (2) with the parameters

$${\mu }_B=0.2999,{\mu }_T=0.8,{\lambda }_B=0.265,{\alpha }_B=0.1,p=0.8,{\beta }_B=0.1,{\beta }_T=0.48,$$

subject to the initial value

$$S_B\left(0\right)=0.1,I_B\left(0\right)=0.2,I_T\left(0\right)=0.1,$$

and$t\in [0,{10}^5]$.

In this case, ${\mathcal{R}}_0=0.2830<1$ and $F_1^{*}=(1,0,0)$ is globally asymptotically stable. We use the Euler scheme, the classical fourth-stages Runge–Kutta (RK4) scheme and NSFD schemes (6) to numerically solve system (2). The numerical solutions obtained by these schemes are presented in Figure 1. From this figure, we see that the numerical solutions obtained by the Euler scheme and the RK4 schemes are not positive, their boundedness is destroyed. The Euler scheme gives the numerical solutions oscillating near the equilibrium points. Meanwhile, the numerical solution obtained by NSFD schemes preserve the essential properties of the model (2). Besides, the NSFD schemes is easily realized. It is the advantages of the NSFD schemes compared with standard difference schemes. Other examples give the same results. This fact completely agrees with several results in the previous works [11,13,14,15,20,21,22,23].

Example 2 (Dynamics of NSFD schemes in the case ${\mathcal{R}}_0\le 1$).

Consider system (2) with the parameters (see in [3])

$${\mu }_B=0.0002999,{\mu }_T=0.001609,{\lambda }_B=0.0265,{\alpha }_B=0.001,p=0.1,{\beta }_B=0.003,{\beta }_T=0.00048.$$

In this case, ${\mathcal{R}}_0=0.33$ and the number ${\phi }^{*}$ in Theorem 4 is ${\phi }^{*}=37.7358$. Therefore, we take

$$\phi \left(h\right)={\displaystyle \frac{1-e^{-0.03h}}{0.03}}.$$

Numerical solutions obtained from this NSFD scheme (6) is depicted in Figure 2. Clearly, the properties of the continuous-time model are preserved.

Example 3 (Dynamics of NSFD schemes in the case ${\mathcal{R}}_0>1$).

Consider system (2) with the parameters (see in [3])

$${\mu }_B=0.0002999,{\mu }_T=0.001609,{\lambda }_B=0.000265,{\alpha }_B=0.001,p=0.1,{\beta }_B=0.006,{\beta }_T=0.00048.$$

In this case, ${\mathcal{R}}_0=67.54$ and the number ${\phi }^{*}$ in Theorem 4 is ${\phi }^{*}=166.67$. Therefore, we take

$$\phi \left(h\right)={\displaystyle \frac{1-e^{-0.005h}}{0.005}}.$$

Numerical solutions obtained from this NSFD scheme (6) is depicted in Figure 3. Similar to Example 2, the properties of the continuous-time model are preserved.

5. Conclusions

As the main conclusion of this study, we obtain the dynamic consistency of NSFD schemes in relation to the continuous-time model. That is, NSFD schemes preserve essential properties of the continuous-time model, while standard finite difference schemes (SFDS) cannot preserve these properties.

On the other hand, it is worth to note that our results constitute a generalization of those in [10].

In future works, the results in this paper will be useful in order to develop other essential applied models, specially those models possessing non-hyperbolic equilibrium points with global stability property.