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

: In this paper, a discrete-time model for Babesiosis disease, given by means of nonstandard ﬁnite difference (NSFD) schemes, is ﬁrst 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.


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 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 R 0 ≤ 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.

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, where the total population of bovine N B (t) is divided into three subpopulations: bovines which may become infected (Susceptible S B (t)), bovines infected by the Babesia parasite (Infected I B (t)), and bovines which have been treated for the Babesiosis (Controlled C B (t)). That is, as in [3], The total population of ticks N T (t) is divided into two subpopulations: ticks which may become infected S T (t) and ticks infected by the Babesia parasite I T (t). As in [3], the parameter µ B represents the birth rate of the bovines and it is assumed equal to their natural death; µ 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 β B ; any susceptible tick can become infected when it stings an infected bovine, at a rate β T ; p is the probability for a susceptible tick to be born from an infected one; λ B represents the fraction of the infected bovines which are controlled, that is, treated against Babesia, while α 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 (2) and the region is a positive invariant set for the system (2). For this system, the threshold parameter is given by Concerning system (2), the following results are obtained.
(ii) If R 0 ≤ 1, then the disease-free point F * 1 is globally asymptotically stable; otherwise, the diseasefree point F * 1 is unstable. (iii) If 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. (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 = ∆t = T/N the time step size of the discretization:

Nonstandard Finite Difference Schemes for System
and let t k = kh for k = 0, 1, . . . , N. We now denote by S k B , I k B and I k T the approximated values of S B (t k ), I B (t k ) and I T (t k ), respectively. Following the Mickens' methodology, we propose NSFD schemes of the form where The conditions on ϕ(h) will be determined so that the properties of system (2) are preserved. For convenience, h is omitted in some presentations.
Proof. The theorem is proved by induction. First, rewrite system (6) in the form Clearly, if (S k B , I k B , I k T ) ∈ Ω and condition (7) Adding the first and second equations of (8), we obtain Thus, the proof is complete.
Similar to Proposition 1 in [3], we obtain the following result. (i) If R 0 ≤ 1, then the disease-free equilibrium point F * 1 is globally asymptotically stable. (ii) If 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 : Ω → R + defined by Clearly, V is continuous, V S k B , I k B , I k T ≥ 0 for all S k B , I k B , I k T ∈ Ω, and V(F * 1 ) = 0. From (6), we have which implies that ∆V ≤ 0 for all S k B , I k B , I k T ∈ Ω. Let G * be the largest positively invariant set contained in Then, by using (9) we have that Consequently, it is easy to verify that F * 1 is G * -globally asymptotically stable if R 0 ≤ 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 R 0 ≤ 1.
Thus, the theorem is proved.
Next, we shall investigate the local stability of system (6) in the case R 0 > 1. Denote by σ(J) the set of eigenvalues of the Jacobian matrix J of system (2) at F * 2 . Based on mathematical analyses in [3], we have if 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., Now we denote by K the Jacobian matrix of system (6) at F * 2 . Then K = I + ϕJ (see [13,14]), where I is 3 × 3 unit matrix. Therefore, the eigenvalues Λ i of J correspond to the eigenvalues Ψ i of K determined by 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 This is equivalent to The above inequality is satisfied if Conversely, if R 0 ≤ 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 R 0 > 1. Suppose Λ 1 , Λ 2 , Λ 3 are the eigenvalues of the Jacobian matrix of system (2) at F * 2 . Set Then, F * 2 is locally asymptotically stable if Summing up the results of this section, we obtain the following theorem.

Remark 1.
Concerning Theorem 2, we have to underline three important observations: • Part (i) of Theorem 2 is only appropriate when Ω is a positively invariant set of (6).

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.
In this case, 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 R 0 ≤ 1). Consider system (2) with the parameters (see in [3]) µ B = 0.0002999, µ T = 0.001609, λ B = 0.0265, α B = 0.001, p = 0.1, β B = 0.003, β T = 0.00048.
In this case, R 0 = 67.54 and the number ϕ * in Theorem 4 is ϕ * = 166.67. Therefore, we take Numerical solutions obtained from this NSFD scheme (6)

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.

Conflicts of Interest:
The authors declare no conflict of interest.

Data Availability:
The data used to support the findings of this study are included within the article.