A dynamically-consistent nonstandard finite difference scheme for the SICA model

In this work, we derive a nonstandard finite difference scheme for the SICA (Susceptible-Infected-Chronic-AIDS) model and analyze the dynamical properties of the discretized system. We prove that the discretized model is dynamically consistent with the continuous, maintaining the essential properties of the standard SICA model, namely, the positivity and boundedness of the solutions, equilibrium points, and their local and global stability.


Introduction
The Human Immunodeficiency Virus (HIV) is responsible for a very high number of deaths worldwide. Acquired ImmunoDeficiency Syndrome (AIDS) is a disease of the human immune system caused by infection with HIV. The HIV virus can be transmitted by several ways but there is no cure or vaccine for AIDS. Nevertheless, antiretroviral (ART) treatment improves health, prolongs life and reduces the risk of HIV transmission. The ART treatment increases life expectation but has some limitations. For instance, it doesn't restore health, has some side effects, and is very expensive. Individuals infected with HIV are more likely to develop tuberculosis because of their immunodeficiency, so a model that considers HIV and tuberculosis is very interesting to investigate. One can find many studies in the literature [1][2][3]. A TB-HIV/AIDS co-infection model, that contains the celebrated SICA (Susceptible-Infected-Chronic-AIDS) model as a sub-model, was first proposed in [4]: I(t) = λ(t)S (t) − (ρ + φ + µ)I(t) + αA(t) + ωC(t), the total population at time t. The meaning of the parameters that appear in the SICA model (1.1)-(1.3) are given in Table 1. The model considers a varying population size in a homogeneously mixing population, subdividing the human population into four mutually-exclusive compartments: -susceptible individuals (S ); -HIV-infected individuals with no clinical symptoms of AIDS (the virus is living or developing in the individuals but without producing symptoms or only mild ones) but able to transmit HIV to other individuals (I); -HIV-infected individuals under ART treatment (the so called chronic stage) with a viral load remaining low (C); -HIV-infected individuals with AIDS clinical symptoms (A).
The SICA model has some assumptions. It can be seen in [5] that the susceptible population is increased by the recruitment of individuals into the population, assumed susceptible at a rate Λ. All individuals suffer from natural death at a constant rate µ. Susceptible individuals S acquire HIV infection, following effective contact with people infected with HIV, at rate λ (1.2), where β is the effective contact rate for HIV transmission. The modification parameter η A ≥ 1 accounts for the relative infectiousness of individuals with AIDS symptoms, in comparison to those infected with HIV with no AIDS symptoms. Individuals with AIDS symptoms are more infectious than HIV-infected individuals because they have a higher viral load and there is a positive correlation between viral load and infectiousness. On the other hand, η C ≤ 1 translates the partial restoration of the immune function of individuals with HIV infection that use correctly ART. The SICA mathematical model (1.1)-(1.3) is well-studied in the literature [5]. It has shown to provide a proper description with respect to the HIV/AIDS situation in Cape Verde [6] and recent extensions include stochastic transmission [7] and fractional versions with memory and general incidence rates [8]. Here our main aim is to propose, for the first time in the literature, a discrete-time SICA model.
For most nonlinear continuous models in engineering and natural sciences, it is not possible to obtain an exact solution [9], so a variety of methods have been constructed to compute numerical solutions [10,11]. It is well known that numerical methods, like the Euler and Runge-Kutta, among others, often fail to solve nonlinear systems. One of the reasons is that they generate oscillations and unsteady states if the time step size decreases to a critical size [12]. Among available approaches to address the problem, the nonstandard finite discrete difference (NSFD) schemes, introduced by Mickens in [13,14], have been successfully applied to several different epidemiological models [15,16]. Precisely, the NSFD schemes were created to eliminate or reduce the occurrence of numerical instabilities that generally arise while using other methods. This is possible because there are some designed laws that systems must satisfy in order to preserve the qualitative properties of the continuous model, such as positivity, boundedness, stability of the equilibrium points, conservation laws, and others [17]. The literature on Mickens-type NSFD schemes is now vast [18,19]. The paper [20] considers the NSFD method of Mickens and apply it to a dynamical system that models the Ebola virus disease. In [21], a NSFD scheme is designed in which the Metzler matrix structure of the continuous model is carefully incorporated and both Mickens' rules on the denominator of the discrete derivative and the nonlocal approximation of nonlinear terms are used. In that work the general analysis is detailed for a MSEIR model. In [22], the authors summarize NSFD methods and compare their performance for various step-sizes when applied to a specific two-sex (male/female) epidemic model; while in [23] it is shown that Mickens' approach is qualitatively superior to the standard approach in constructing numerical methods with respect to productive-destructive systems (PDS's). NSFD schemes for PDS's are also investigated in [24]; NSFD methods for predator-prey models with the Beddington-De Angelis functional response are studied in [25]. Here we propose and investigate, for the first time in the literature, the dynamics of a discretized SICA model using the Mickens NSFD scheme.
The paper is organized as follows. Some considerations, regarding the continuous SICA model and the stability of discrete-time systems, are presented in section 2. The original results are then given in section 3: we start by introducing the discretized SICA model; we find the equilibrium points, prove the positivity, Theorem 3, and boundedness of the solutions, Theorem 4; we also establish the local stability of the disease free equilibrium point of the discrete model, Theorem 5, as well as the global stability of the equilibrium points, Theorems 6 and 7. In section 4, we provide some numerical simulations to illustrate the stability of the NSFD discrete SICA model using a case study. We end with section 5 of conclusion.

Preliminaries
In this section, we collect some preliminary results about the continuous SICA model [4], as well as results for the stability of discrete-time systems [26], useful in our work.

The continuous SICA model
All the information in this section is proved in [4]. Each solution (S (t), I(t), C(t), A(t)) of the continuous model much satisfy S (0) ≥ 0, I(0) ≥ 0, C(0) ≥ 0, and A(0) ≥ 0, because each equation represents groups of human beings. Adding the four equations of (1.1), one has so that Therefore, the biologically feasible region is given by which is positively invariant and compact. This means that it is sufficient to study the qualitative dynamics in Ω. The model has two equilibrium points: a disease free and an endemic one. The disease free equilibrium (DFE) point is given by Following the approach of the next generation matrix [27], the basic reproduction number R 0 for model (1.1), which represents the expected average number of new infections produced by a single HIVinfected individual in contact with a completely susceptible population, is given by where along all the manuscript we use C 1 = ρ + φ + µ, C 2 = α + µ + d, and C 3 = ω + µ. The endemic point has the following expression: where D = −(λ * + µ)(µ(C 3 (ρ + C 2 ) + C 2 φ + ρd) + ρωd) and which is positive if R 0 > 1. The explicit expression of the endemic equilibrium point of (1.1) is given by Regarding the stability of the equilibrium points, Theorem 3.1 and Proposition 3.4 of [4] establish the persistence of the endemic point. The disease is persistent in the population if the infected cases with AIDS are bounded away from zero or the population S disappears. The local stability of the endemic point is given in Theorem 3.8 of [4]. Lemma 3.5 of [4] states that the DFE is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1. Finally, Theorem 3.6 of [4] asserts that, under suitable conditions, the DFE point is globally asymptotically stable. Here we prove similar properties in the discrete-time setting (section 3). For that we now recall an important tool for difference equations.

The Schur-Cohn criterion
One of the main tools that provides necessary and sufficient conditions for the zeros of a nth-degree polynomial to lie inside the unit disk is the Schur-Cohn criterion [26]. This result is useful for studying the stability of the zero solution of a kth-order difference equation or to investigate the stability of a k-dimensional system of the form where p(λ) in (2.3) is the characteristic polynomial of the matrix A. Let us introduce some preliminaries before presenting the Schur-Cohn criterion. Namely, let us define the inners of a matrix B = (b i j ). The inners of a matrix are the matrix itself and all the matrices obtained by omitting successively the first and last columns and first and last rows. A matrix B is said to be positive innerwise if the determinant of all its inners are positive.
Theorem 1 (The Shur-Cohn criterion [26]). The zeros of the characteristic polynomial (2.3) lie inside the unit disk if, and only if, the following holds: Using the Schur-Cohn criterion, one may obtain necessary and sufficient conditions on the p i coefficients such that the zero solution of (2.3) is locally asymptotically stable.

Main results
We begin by proposing a discrete-time SICA model.

The NSFD scheme
One of the important features of the discrete-time epidemic models obtained by Mickens method is that they present the same features as the corresponding original continuous-time models. Here, we construct a dynamically consistent numerical NSFD scheme for solving (1.1) based on [13,14,17]. Let us define the time instants t n = nh with n integer, h = t n+1 − t n as the time step size, and (S n , I n , C n , A n ) as the approximated values of (S (nh), I(nh), C(nh), A(nh)). Thus, the NSFD scheme for model (1.1) takes the following form: The nonstandard schemes are based in two fundamental principles [14,17]: 1. Regarding the first derivative, we have where ν(h) and ψ(h) are the numerator and denominator functions that satisfy the requirements In general, the numerator function can be selected to be ν(h) = 1. We will make this choice here. Generally, the denominator function is nontrivial. Based on Mickens work [28,29], when we write explicitly S n+1 using φ(h) = h, we have in the denominator the term 1 + µh, which means that we can use ψ(h) = e µh −1 µ as the denominator function. Throughout our study, ψ(h) = e µh −1 µ but, for brevity, we write ψ(h) = ψ. 2. Both linear and nonlinear functions of x(t) and its derivatives may require a "nonlocal" discretization. For example, x 2 can be replaced by x k x k+1 .
Since model (3.1) is linear in S n+1 , I n+1 , C n+1 , and A n+1 , we can obtain their explicit form:

Positivity of solutions
The first result to be shown is the positivity of the solutions. Proof. Let us assume that S (0), I(0), C(0), A(0) are positive. We only need to show that I n+1 is positive. The denominator is given by which can be rewritten as and, simplifying, we get Since all parameter values are positive, (S n+1 , I n+1 , C n+1 , A n+1 ) is positive for all n ≥ 0.

Conservation law
The second result to be shown is the conservation law or boundedness of the solutions.
Theorem 4. The NSFD scheme defines the discrete dynamical system (3.2) oñ Ω = (S n , I n , C n , A n ) : 0 ≤ S n + I n + C n + A n ≤ Λ µ . (3.4) Proof. Let the total population be N n = S n + I n + C n + A n . Adding the four equations of (3.1), we have By the discrete Grownwall inequality, if 0 < N(0) < Λ µ , then and, since 1 1 + µψ < 1, we have N n → Λ µ as n → ∞. We conclude that the feasible regionΩ is maintained within the discrete scheme.

Elementary stability
A difference scheme that approximates a first-order differential system is elementary stable if, for any value of the step size, its fixed-points are exactly those of the differential system. Furthermore, when applied to the associated linearized differential system, the resulting difference scheme has the same stability/instability properties [30].
The continuous and discrete system have the same equilibria. The disease free equilibrium (DFE) point is given by The explicit expression of the endemic equilibrium point of (3.2) can be rewritten as the one of the continuous case: when we obtain the endemic equilibrium point. The explicit expression of the endemic equilibrium point of (3.2) can be rewritten as the one of the continuous case: After some computations, we get the following equalities:

Local stability of the DFE point
Let us discuss the stability of the proposed NSFD scheme at the DFE point E 0 (3.6).
Remark 1. Several articles use the next-generation matrix approach presented in [31]. For that, however, the matrices F + T must be non-negative. Our model does not satisfy such condition.
The following technical lemma has an important role in our proofs.
Proof. The linearization of (3.2) at the DFE E 0 is: The characteristic polynomial of J(E 0 ) has the following expression: where p 3 (λ) is given by that is, Since p 4 (1) > 0, we have (−µ − 1)p 3 (1) > 0, and recalling that µ > 0, we conclude that Therefore, We conclude that the first condition of Theorem 1 is satisfied if R 0 < 1.
, p 2 < 1 + p 4 , and then the third condition of Theorem 1 is satisfied.
We are now in condition to prove the main result of this section.
(1−C 2 )(1−C 3 ) , p 2 < 1 + p 4 , and (3.9) and (3.12) are satisfied, then, provided R 0 < 1, the disease free equilibrium point of the discrete system (3.2) is locally asymptotically stable. If the previous conditions are not satisfied, then the disease free equilibrium point is unstable.
Proof. The result follows by Theorem 1 and Propositions 1, 2 and 3. If any of the conditions enumerated are not satisfied, at least one of the roots of the characteristic polynomial lies outside the unit circle, so the disease free equilibrium point is unstable.

Global stability of the equilibrium points
Now we prove that R 0 is a critical value for global stability: when R 0 < 1, the disease free equilibrium point is globally asymptotically stable; when R 0 > 1, the endemic equilibrium is globally asymptotically stable.
Theorem 6. If R 0 < 1, then the DFE point of the discrete-time SICA model (3.1) is globally asymptotically stable.
Theorem 7. If R 0 > 1, then the endemic equilibrium point of the discrete-time SICA model (3.1) is globally asymptotically stable.
Proof. We construct a sequence {Ṽ(n)} +∞ n=1 of the form where g(x) = x − 1 − ln(x), x ∈ R + . Clearly, g(x) ≥ 0 with equality holding true only if x = 1. We have Similarly, The difference ofṼ(n) satisfies Therefore, {Ṽ(n)} +∞ n=1 is a monotone decreasing sequence for any n ≥ 0. SinceṼ(n) ≥ 0 and lim n→∞ Ṽ (n + 1) −Ṽ(n) = 0, we obtain that lim n→∞ S n+1 = S * , lim n→∞ I n+1 = I * , lim n→∞ C n+1 = C * and lim n→∞ infected individuals under ART treatment have a very low probability of transmitting HIV [40]. For the parameter η A ≥ 1, which accounts the relative infectiousness of individuals with AIDS symptoms, in comparison to those infected with HIV with no AIDS symptoms, we assume, based on [41], that η A = 1.35. For (η C , η A ) = (0.04, 1.35), the estimated value of the HIV transmission rate is equal to β = 0.695. Using these parameter values, the basic reproduction number is R 0 = 4.5304 and the endemic equilibrium point (S * , I * , C * , A * ) = (145276, 48136.4, 461146, 3580.57). In Figure 1, we show graphically the cumulative cases of infection by HIV/AIDS in Cape Verde given in Table 2, together with the curves obtained from the continuous-time SICA model (1.1) and our discrete-time SICA model (3.1). Our simulations of the continuous and discrete models were done with the help of the software Wolfram Mathematica, version 12.1. For the continuous model, we have used the command NSolve, that computes the solution by interpolation functions. Our implementation for the discrete case makes use of the Mathematica command RecurrenceTable. To illustrate the global stability of the endemic equilibrium (EE), predicted by our Theorem 7, we consider different initial conditions borrowed from [5], from different regions of the plane:  The obtained results are given in Figure 2.

Conclusions
In this work, we proposed a discrete-time SICA model, using Mickens' nonstandard finite difference (NSFD) scheme. The elementary stability was studied and the global stability of the equilibrium points proved. Finally, we made some numerical simulations, comparing our discrete model with the continuous one. For that, we have used the same data, following the case study of Cape Verde. Our conclusion is that the discrete model can be used with success to describe the reality of Cape Verde, as well as to properly approximate the continuous model. All our simulations have been done using the numerical computing environment Mathematica, version 12.1, running on an Apple MacBook Pro i5 2.5 GHz with 16Gb of RAM. The solutions of the models were found in "real time".
Mickens was a pioneer in NSFD schemes. Throughout the years, other NSFD schemes were developed. Roughly speaking, different Mickens-type methods differ on the denominator functions and the discretization, depending on concrete conditions that the continuous model under study must satisfy. In [20], the incidence rate is combined, while in [21] all parameters are constant. Other types of NSFD are presented, e.g., in [22][23][24][25], which can be used if the system satisfy some conditions and a suitable denominator function is constructed. For such schemes, the incidence functions are different from ours. In [23], for example, it is fundamental to rewrite the system and the discretization method and the denominator function are different from the ones we use here. The article [24] uses the same approach of [23] and the model has a bilinear incidence function. Positive and elementary stable nonstandard finite-difference methods are also considered in [25]. Mickens set the field, but several different authors developed and are developing other related discretization methods. For the SICA model, however, as we have shown, a new NSFD scheme is not necessary and the standard Mickens' method provides a well posed discrete-time model with excellent results, without the need to impose additional conditions to the model.
Interior (CMA-UBI), project UIDB/00212/2020; Delfim F. M. Torres through the Center for Research and Development in Mathematics and Applications (CIDMA), project UIDB/04106/2020. They are very grateful to three anonymous referees, who kindly reviewed an earlier version of this manuscript and provided valuable suggestions and comments.