**1. Introduction**

History has recorded that infectious diseases have caused devastation in the human population. Despite the great advances in epidemic control, it was believed that infectious diseases would soon be eradicated, but this has clearly not been the case. Microorganisms adapt and evolve, and consequently, new infectious diseases such as AIDS, Ebola or COVID-19 appear, which cause many deaths. In addition, the genome of some microorganisms can sometimes change slightly and consequently, they can acquire resistance to some drugs [1]. According to the World Health Organization, since its first registration in the 1980s, Human Immunodeficiency Virus (HIV), the causative agent of Acquired Immunodeficiency Syndrome (AIDS), has caused more than 35 million deaths worldwide [2]. The greatest impact of deaths caused by AIDS occurs in underdeveloped or very poor countries, especially in sub-Saharan Africa [2,3].

HIV is an RNA virus that belongs to the retroviridae family, specifically to the lentivirus subfamily, and acts against the immune system, weakening its defense systems against infections and certain types of cancer, which is why the infected person gradually loses its immunodeficiency [4,5]. The HIV replication process is active and dynamic in the sense that when the virus enters the body, the cells that have the CD4+ receptor are infected, most of them are TCD4+ lymphocytes. After entering the cell, the HIV virus can remain latent, replicate in a controlled manner, or undergo massive replication that results in a cytopathic effect for the infected cell. In most lymphocytes the virus is latent, and the infection gradually decreases the amount of these in both the tissues and the blood. This leads the patient to a severe state of cellular immunosuppression, and then a group of microorganisms causes infections. As a consequence, there is a great mortality of people affected by HIV [6].

**Citation:** Arenas, A.J.; González-Parra, G.; Naranjo, J.J.;

Cogollo, M.; De La Espriella, N. Mathematical Analysis and Numerical Solution of a Model of HIV with a Discrete Time Delay. *Mathematics* **2021**, *9*, 257. http://doi.org/10.3390/math9030257

Academic Editor: Francisco Rodríguez Received: 31 December 2020 Accepted: 25 January 2021 Published: 28 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Epidemiologists conduct scientific experiments, sometimes in controlled settings through self-experimentation, to analyze the spread and possible control strategies of infectious diseases. However, designing such controlled experiments is sometimes impossible due to ethical concerns and the possible collection of erroneous data [1,7,8]. These reasons motivate the possibility of using mathematical models as tools to corroborate the perception of disease transmission, test theories, and suggest better intervention and control strategies.

Recently, there have been a growing literature regarding mathematical models for virus infection within-host [9–19]. These mathematical models include a variety of characteristics related to the viral dynamics. For instance, some models include discrete time delays, cell to cell viral transmissions, and the most well-known virus to cell transmission [10,15,17,19–21]. This article presents a new mathematical model, by means of a set of differential equations with delay, to determine the effect of how to produce viruses by target cells inside the dynamics of viruses. In this case, two types of virus-infected cells are analyzed: the cells in the eclipse phase that are not producing the virus *IE*, and the cells that are actively producing the virus *I*. The cells in the eclipse phase change to the state of virus production at a *m* rate, and the mortality rate of each cell type is *δIE* and *δI*, respectively. Cells in the eclipse phase may die because they could be recognized as infectious by mediators of innate immunity or due to the accumulation of DNA intermediates in the cell cytoplasm, see [22,23].

Numerous mathematical models represented by means of a system of differential equations, with or without delay, have been discretized by means of the non-standard finite difference method proposed by Ronald Mickens, see [24–34]. Their use is mainly because they are very effective in preserving certain qualitative properties of the original differential equations and the convergence, consistency and stability of their solutions have been demonstrated, see [35–46].

In this study, we also design a non-standard finite difference (*NSFD*) scheme that allows us to obtain numerical solutions of a set of delayed and ordinary differential equations, which describes the dynamics of HIV infection within-host. First, we apply the techniques designed by Mickens for the construction of the non-standard finite difference (*NSFD*) scheme to our HIV mathematical model. Secondly, we prove some main properties of the non-standard finite difference (*NSFD*) scheme, and that agree with qualitative properties on the HIV mathematical model with discrete time delay. One important property that the non-standard finite difference (*NSFD*) scheme has is that it allows us to guarantee accurate and positive solutions. This is very important when solving inverse problems related to estimation of parameters [13,47–51]. Finally, we perform some numerical simulations that show the advantages regarding accuracy and computational cost.

This paper is organized as follows. In Section 2, we present the mathematical model of HIV within-host with discrete time delay. Section 3 is devoted to the stability mathematical analysis. In Section 4, we construct the *NSFD* numerical scheme. We include in this section the study of some properties of this numerical scheme such as stability analysis. In Section 5, the numerical simulation results using the constructed *NSFD* scheme are shown, and the last section is devoted to the conclusions.

#### **2. Mathematical Model of HIV Within-Host with Discrete Time Delay**

We construct the mathematical model using a combination of virus facts, hypotheses, and previous proposed mathematical models within-host [9,11,14,52–55]. Despite there has been a growing number of studies related to viral dynamics within-host there are still some aspects that are not well understood [8,9,18,56–58]. Moreover, mathematical models include assumptions that make them more tractable to be able to extract useful information and test different hypotheses [9,55,59,60]. We start noticing that it has been argued that at least in vitro, most HIV-infected cells die before virus production begins [22,61,62]. Virusproducing cells produce virions *V* at a rate of *NδI*, where *N* is the average number of infectious virions released by an infected cell during its lifetime. In general, it is accepted that most virions produced by infected cells are not infectious [63–65], and since these

virions are not contributing to the infection of new cells, non-infectious virions are not considered in this model. On the other hand, *V* virions that are infectious can be removed by the immune system from the population of virus-free cells at intrinsic clearance rate *C*, or they can infect target cells (CD4+ T) at a *β* rate, with *T* the concentration of target cells, where Λ is the generation rate of uninfected CD4+ T cells and *μ*<sup>0</sup> mortality rate of uninfected cells. In the constructed mathematical model there are two classes of infected cells. The first one is the class that includes the cells in the eclipse phase which are not making the virus, *IE*. The second class include the cells that are actively producing the virus, *I*. Cells in the eclipse phase transition to the class *I* at a rate, *m*. These cells die at rate *δIE* . The cells in class *I* then die at rate *δI*. Cells in the eclipse phase die due to the immune systems. Notice that the number of targets cells *T* are not virus-infected and vary depending on the parameters Λ and the particular death rate for target cells *μ*0.

Based on the previous assumptions, we propose a model that describes the intracellular dynamics of HIV and is given by the following system of ordinary differential equations,

$$\begin{aligned} \frac{dI\_E(t)}{dt} &= \beta T(t)V(t) - (m + \delta\_{I\_E})I\_E(t) \\ \frac{dI(t)}{dt} &= mI\_E(t) - \delta\_I I(t) \\ \frac{dV(t)}{dt} &= N\delta\_I I(t) - CV(t) - \beta T(t)V(t) \\ \frac{dT(t)}{dt} &= \Lambda - \beta T(t)V(t) - \mu\_0 T(t) \end{aligned} \tag{1}$$

Notice, that the transmission term *βT*(*t*)*V*(*t*) also appears in the equation for virions because of the assumption that it takes only one virion to infect a target cell [52,66]. The possibility of multiple infected cells is excluded [66]. It also should be noted that during the eclipse phase (the time from viral entry to the active production of viral particles) the infected cells are not producing virions [8,9,16,67,68]. This delay affects the maximum viral load time and the probability that a viral infection will be established [9,68–70], and therefore should be explicitly modeled [9,68]. For this case, let Δ > 0 be the duration of the eclipse phase and *iE*(*t*, *τ*) the density of cells at a time *t* that were infected *τ* units of time before of *t*, i.e., are infected cells of age *τ*. Then *iE*(*t*, Δ) represents the proportion of cells in the eclipse phase that go into the state of virus production, whereby

$$\frac{dI(t)}{dt} = i\_E(t, \Delta) - \delta\_I I(t). \tag{2}$$

Since the mortality rate of eclipse cells *δIE* is constant, it is appropriate to assume that *iE*(*t*, *τ*) satisfies the equation of McKendrick-Von Foerster age-structured population dynamics model, see [71] that is

$$\frac{\partial \dot{i}\_E(t, \tau)}{\partial t} + \frac{\partial \dot{i}\_E(t, \tau)}{\partial \tau} = -\delta\_{I\_E} \dot{i}\_E(t, \tau),\tag{3}$$

subject to the boundary condition *iE*(*t*, 0) = *βT*(*t*)*V*(*t*). Thus, the solution of the Equation (3) is given by *iE*(*t*, *τ*) = *βT*(*t* − *τ*)*V*(*t* − *τ*)*e* <sup>−</sup>*δIE <sup>τ</sup>*. Therefore, Equation (2) is given by

$$\frac{dI(t)}{dt} = \beta T(t-\Delta)V(t-\Delta)e^{-\delta\_{l\_{\mathbb{E}}}\Delta} - \delta\_{I}I(t). \tag{4}$$

The total number of cells in the eclipse phase is given by *IE*(*t*) = - Δ <sup>0</sup> *iE*(*t*, *τ*)*dτ*, and by integrating both sides of Equation (3) on the interval [0, Δ] we have

$$\frac{dI\_E(t)}{dt} = \beta T(t)V(t) - \beta T(t-\Delta)V(t-\Delta)e^{-\delta\_{I\_E}\Delta} - \delta\_{I\_E}I\_E(t). \tag{5}$$

Thus the mathematical model given in (1) takes the form

$$\begin{aligned} \frac{dI\_E(t)}{dt} &= \beta T(t)V(t) - \beta T(t-\Delta)V(t-\Delta)e^{-\delta\_{I\_E}\Delta} - \delta\_{I\_E}I\_E(t) \\ \frac{dI(t)}{dt} &= \beta T(t-\Delta)V(t-\Delta)e^{-\delta\_{I\_E}\Delta} - \delta\_I I(t) \\ \frac{dV(t)}{dt} &= N\delta\_I I(t) - CV(t) - \beta T(t)V(t) \\ \frac{dT(t)}{dt} &= \Lambda - \beta T(t)V(t) - \mu\_0 T(t), \end{aligned} \tag{6}$$

where Δ is the duration of the eclipse phase, *e* <sup>−</sup>*δIE*<sup>Δ</sup> represents the probability that an infected cell will survive a time Δ after viral entry. Notice that Δ is a fixed time delay, and then we have a differential equation with a discrete time delay. For a better understanding of the mathematical model, we can analyze it from the following flow chart shown in Figure 1 [54]:

**Figure 1.** Flow chart of model (6).

In addition, the system (6) satisfies the initial conditions given by

$$T(\mathbf{s}) = \mathfrak{J}\_1(\mathbf{s}), V(\mathbf{s}) = \mathfrak{J}\_2(\mathbf{s}), \mathbf{s} \in [-\Delta, \mathbf{0}], \tag{7}$$

with *<sup>ξ</sup>*1(*s*), *<sup>ξ</sup>*2(*s*) positive continuous functions defined from the interval [−Δ, 0] to <sup>R</sup><sup>2</sup> +, and equipped with the norm *ξ*1,2 = sup −Δ≤*s*≤0 |*ξ*1,2|.

## **3. Properties of the Solutions of the Mathematical Model**

Using the fundamental theory of functional differential equations [72,73], it follows that the solution of the system (6) with the initial condition (7) exists for all *t* ≤ 0 and is unique.

Next, we will establish different dynamic properties of the solution of the mathematical model described by the system of Equation (6). Since the system (6) represents a biological model, it is important to determine the nature of the solution. Thus, if it is assumed that all the parameters are non-negative as well as the initial conditions *IE*(*s*), *I*(*s*), *V*(*s*), *T*(*s*), for *s* ∈ [−Δ, 0]. Therefore, we must guarantee the positivity and boundedness of the solution (*IE*(*t*), *I*(*t*), *V*(*t*), *T*(*t*)) of the system (6) at [0, ∞). The following results characterizes these properties.

**Theorem 1.** *If the initial conditions IE*(0) = *IE*<sup>0</sup> , *I*(0) = *I*0, *V*(0) = *V*0, *T*(0) = *T*<sup>0</sup> *of the mathematical model (6) are positive, then the solutions* (*IE*(*t*), *I*(*t*), *V*(*t*), *T*(*t*)) *of the system (6) are positive for all t* ∈ [0, ∞).

**Proof.** Let us start by noting that the solutions of the differential equations of *IE*(*t*) and *I*(*t*) given in (6) can be written as

$$I\_E(t) = e^{-\delta\_{I\_E}t} \left[ I\_{E\_0} + \int\_{t-\Lambda}^t \beta T(s) V(s) e^{\delta\_{I\_E}s} ds \right],\tag{8}$$

$$I(t) = e^{-\delta\_l t} \left[ I\_0 + \int\_0^t \beta T(s - \Delta) V(s - \Delta) e^{(\delta\_l - \delta\_{l\_E})s} ds \right]. \tag{9}$$

Therefore, the positivity of the solutions *T*(*t*) and *V*(*t*) for all *t* > 0, allows us to guarantee the positivity of *IE*(*t*) and *I*(*t*) and thus of system (6). Thus, for *T*(*t*) given as in (6) we have that *T*(*t*) > 0, for all *t* ≥ 0. Indeed, suppose that the positivity does not holds, therefore there must be a *<sup>t</sup>*<sup>0</sup> <sup>&</sup>gt; 0 such that *<sup>T</sup>*(*t*0) = 0, *dT*(*t*0) *dt* <sup>≤</sup> 0 and *<sup>T</sup>*(*t*) <sup>&</sup>gt; 0 for all *t* ∈ [0, *t*0), because the initial condition *T*<sup>0</sup> > 0. Thus, *T*(*t*) must be negative from some *t*0. However, in the interval [0, *t*0) the function *T*(*t*) is positive, and at point *t*<sup>0</sup> the derivative at *t*<sup>0</sup> is non-positive. Thus, from the fourth equation of model (6), it follows that for *t*0,

$$\frac{dT(t\_0)}{dt} = \Lambda - \beta T(t\_0)V(t\_0) - \mu\_0 T(t\_0) = \Lambda > 0,$$

which contradicts that *dT*(*t*0) *dt* <sup>≤</sup> 0. Therefore, we must have *<sup>T</sup>*(*t*) <sup>&</sup>gt; 0, for all *<sup>t</sup>* <sup>≥</sup> 0. Now, we affirm that if *V*(*t*) is given by the system (6), it follows that

$$V(t) > 0,\text{ for all } t \ge 0. \tag{10}$$

Indeed, suppose that there exists a *<sup>t</sup>*<sup>1</sup> <sup>&</sup>gt; 0 such that *<sup>V</sup>*(*t*1) = 0, *dV*(*t*1) *dt* <sup>≤</sup> 0 and *V*(*t*) > 0 for all *t* ∈ [0, *t*1). Then it holds that *I*(*t*) > 0 for all *t* ∈ [0, *t*1). Otherwise, there should be a *<sup>t</sup>*<sup>2</sup> such that 0 <sup>&</sup>lt; *<sup>t</sup>*<sup>2</sup> <sup>&</sup>lt; *<sup>t</sup>*1, *<sup>I</sup>*(*t*2) = 0, *dI*(*t*2) *dt* <sup>≤</sup> 0 and *<sup>I</sup>*(*t*) <sup>&</sup>gt; 0 for all *<sup>t</sup>* <sup>∈</sup> [0, *<sup>t</sup>*2). Thus, from the second equation of system (6), if follows that for *t*<sup>2</sup> ∈ (Δ, *t*1) it holds

$$\frac{dI(t\_2)}{dt} = \beta T(t\_2 - \Delta)V(t\_2 - \Delta)e^{-\delta\_{l\_\Sigma}\Delta} - \delta\_I I(t\_2) = \beta T(t\_2 - \Delta)V(t\_2 - \Delta)e^{-\delta\_{l\_\Sigma}\Delta}.$$

But, −Δ < *t*<sup>2</sup> − Δ < *t*<sup>1</sup> − Δ < *t*1. From the initial conditions given by (7) and the hypothesis for *<sup>V</sup>*(0) <sup>&</sup>gt; 0 it follows that *<sup>V</sup>*(*t*<sup>2</sup> <sup>−</sup> <sup>Δ</sup>) <sup>&</sup>gt; 0. This, contradicts that *dI*(*t*2) *dt* <sup>≤</sup> 0. Thus, *I*(*t*) > 0, for all *t* ∈ [0, *t*1). Next, from third equation of system (6) for *t*<sup>1</sup> > 0,

$$\frac{dV(t\_1)}{dt} = N\delta\_I I(t\_1) - CV(t\_1) - \beta T(t\_1)V(t\_1) = N\delta\_I I(t\_1) > 0.$$

This is a contradiction since *dV*(*t*1) *dt* <sup>≤</sup> 0. Therefore, *<sup>V</sup>*(*t*) <sup>&</sup>gt; 0 for all *<sup>t</sup>* <sup>≥</sup> 0.

**Theorem 2.** *The solution* (*IE*(*t*), *I*(*t*), *V*(*t*), *T*(*t*)) *of system (6) is uniformly bounded in* [0, ∞)*.*

**Proof.** From system (6), one can see that

$$\frac{dI\_E(t)}{dt} + \frac{dI(t)}{dt} + \frac{dT(t)}{dt} = \Lambda - \delta\_{I\_E} I\_E(t) - \delta\_I I(t) - \mu\_0 T(t) \le \Lambda - M(I\_E(t) + I(t) + T(t)), \quad t \ge 0$$

where *M* = *min δIE* , *δI*, *μ*<sup>0</sup> . This implies that

$$I\varepsilon(t) + I(t) + T(t) \le e^{-Mt} \left[ I\varepsilon\_0 + I\alpha + T\_0 + \int\_0^t \Lambda e^{Ms} ds \right] \\ = e^{-Mt} \left[ I\varepsilon\_0 + I\alpha + T\_0 \right] + \frac{\Lambda}{M} \left( 1 - e^{-Mt} \right).$$

Thus,

$$\limsup\_{t \to \infty} (I\_E(t) + I(t) + T(t)) \le \frac{\Lambda}{M}.$$

Accordingly, *IE*(*t*), *I*(*t*) and *T*(*t*) are uniformly boundedness. Even more, given *ε* > 0, there exists *t*<sup>1</sup> > 0 such that for all *t* ≥ *t*1,

$$I(t) \le \frac{\Lambda}{M} + \varepsilon\_{-}$$

Then, since *dV*(*t*) *dt* <sup>≤</sup> *<sup>N</sup>δ<sup>I</sup> <sup>I</sup>*(*t*) <sup>−</sup> *CV*(*t*), and *<sup>V</sup>*(*t*) <sup>&</sup>gt; 0 for all *<sup>t</sup>* <sup>&</sup>gt; *<sup>t</sup>*1, one obtains that

$$\frac{dV(t)}{dt} + \mathcal{C}V(t) \le N\delta\_I \left(\frac{\Lambda}{M} + \varepsilon\right).$$

It follows that

$$V(t) \le V\_0 \varepsilon^{-Ct} + \left(1 - \varepsilon^{-Ct}\right) \left(\frac{N\delta\_I}{C}\right) \left(\frac{\Lambda}{M} + \varepsilon\right).$$

As *t* −→ ∞, then *V*(*t*) ≤ *Nδ<sup>I</sup> C*  Λ *<sup>M</sup>* <sup>+</sup> *<sup>ε</sup>* . Since > 0, *V*(*t*) is uniformly boundedness. This completes the proof.
