Mathematical Analysis and Numerical Solution of a Model of HIV with a Discrete Time Delay

Abstract

We propose a mathematical model based on a set of delay differential equations that describe intracellular HIV infection. The model includes three different subpopulations of cells and the HIV virus. The mathematical model is formulated in such a way that takes into account the time between viral entry into a target cell and the production of new virions. We study the local stability of the infection-free and endemic equilibrium states. Moreover, by using a suitable Lyapunov functional and the LaSalle invariant principle, it is proved that if the basic reproduction ratio is less than unity, the infection-free equilibrium is globally asymptotically stable. In addition, we designed a non-standard difference scheme that preserves some relevant properties of the continuous mathematical model.

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].

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,10,11,12,13,14,15,16,17,18,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,20,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 $I_E,$ 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 ${\delta }_{I_E}$ and ${\delta }_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,25,26,27,28,29,30,31,32,33,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,36,37,38,39,40,41,42,43,44,45,46].

In this study, we also design a non-standard finite difference $\left(NSFD\right)$ 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 $\left(NSFD\right)$ scheme to our HIV mathematical model. Secondly, we prove some main properties of the non-standard finite difference $\left(NSFD\right)$ 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 $\left(NSFD\right)$ 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,48,49,50,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,53,54,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,57,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]. Virus-producing cells produce virions V at a rate of $N{\delta }_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,64,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 $\beta$ rate, with T the concentration of target cells, where $\mathsf{\Lambda }$ is the generation rate of uninfected CD4+ T cells and ${\mu }_0$ 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, $I_E$. 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 ${\delta }_{I_E}$. The cells in class I then die at rate ${\delta }_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 $\mathsf{\Lambda }$ and the particular death rate for target cells ${\mu }_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{array}{cc}\hfill \frac{d{I}_E\left(t\right)}{dt}& =\beta T\left(t\right)V\left(t\right)-(m+{\delta }_{I_E})I_E\left(t\right)\hfill \\ \hfill \frac{dI\left(t\right)}{dt}& =m{I}_E\left(t\right)-{\delta }_{I}I\left(t\right)\hfill \\ \hfill \frac{dV\left(t\right)}{dt}& =N{\delta }_{I}I\left(t\right)-CV\left(t\right)-\beta T\left(t\right)V\left(t\right)\hfill \\ \hfill \frac{dT\left(t\right)}{dt}& =\mathsf{\Lambda }-\beta T\left(t\right)V\left(t\right)-{\mu }_{0}T\left(t\right).\hfill \end{array}$$

(1)

Notice, that the transmission term $\beta T\left(t\right)V\left(t\right)$ 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,69,70], and therefore should be explicitly modeled [9,68]. For this case, let $\Delta >0$ be the duration of the eclipse phase and $i_E(t,\tau )$ the density of cells at a time t that were infected $\tau$ units of time before of t, i.e., are infected cells of age $\tau$. Then $i_E(t,\Delta )$ represents the proportion of cells in the eclipse phase that go into the state of virus production, whereby

$$\frac{dI\left(t\right)}{dt}=i_E(t,\Delta )-{\delta }_{I}I\left(t\right).$$

(2)

Since the mortality rate of eclipse cells ${\delta }_{I_E}$ is constant, it is appropriate to assume that $i_E(t,\tau )$ satisfies the equation of McKendrick-Von Foerster age-structured population dynamics model, see [71] that is

$$\frac{\partial i_E(t,\tau )}{\partial t}+\frac{\partial i_E(t,\tau )}{\partial \tau }=-{\delta }_{I_E}{i}_E(t,\tau ),$$

(3)

subject to the boundary condition $i_E(t,0)=\beta T\left(t\right)V\left(t\right)$. Thus, the solution of the Equation (3) is given by $i_E(t,\tau )=\beta T(t-\tau )V(t-\tau )e^{-{\delta }_{I_E}\tau }.$ Therefore, Equation (2) is given by

$$\frac{dI\left(t\right)}{dt}=\beta T(t-\Delta )V(t-\Delta )e^{-{\delta }_{I_E}\Delta }-{\delta }_{I}I\left(t\right).$$

(4)

The total number of cells in the eclipse phase is given by $I_E\left(t\right)={\int }_0^{\Delta }{i}_E(t,\tau )d\tau$, and by integrating both sides of Equation (3) on the interval $[0,\Delta ]$ we have

$$\begin{array}{c}\frac{d{I}_E\left(t\right)}{dt}=\beta T\left(t\right)V\left(t\right)-\beta T(t-\Delta )V(t-\Delta )e^{-{\delta }_{I_E}\Delta }-{\delta }_{I_E}{I}_E\left(t\right).\end{array}$$

(5)

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

$$\begin{array}{cc}\hfill \frac{d{I}_E\left(t\right)}{dt}& =\beta T\left(t\right)V\left(t\right)-\beta T(t-\Delta )V(t-\Delta )e^{-{\delta }_{I_E}\Delta }-{\delta }_{I_E}{I}_E\left(t\right)\hfill \\ \hfill \frac{dI\left(t\right)}{dt}& =\beta T(t-\Delta )V(t-\Delta )e^{-{\delta }_{I_E}\Delta }-{\delta }_{I}I\left(t\right)\hfill \\ \hfill \frac{dV\left(t\right)}{dt}& =N{\delta }_{I}I\left(t\right)-CV\left(t\right)-\beta T\left(t\right)V\left(t\right)\hfill \\ \hfill \frac{dT\left(t\right)}{dt}& =\mathsf{\Lambda }-\beta T\left(t\right)V\left(t\right)-{\mu }_{0}T\left(t\right),\hfill \end{array}$$

(6)

where $\Delta$ is the duration of the eclipse phase, $e^{-{\delta }_{I_E}\Delta }$ represents the probability that an infected cell will survive a time $\Delta$ after viral entry. Notice that $\Delta$ 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]:

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

$$\begin{array}{c}\hfill T\left(s\right)={\xi }_1\left(s\right),V\left(s\right)={\xi }_2\left(s\right),s\in [-\Delta ,0],\end{array}$$

(7)

with ${\xi }_1\left(s\right),{\xi }_2\left(s\right)$ positive continuous functions defined from the interval $[-\Delta ,0]$ to ${\mathbb{R}}_{+}^2$, and equipped with the norm $\parallel {\xi }_{1,2}\parallel =\underset{-\Delta \le s\le 0}{sup}\left|{\xi }_{1,2}\right|.$

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\le 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 $I_E\left(s\right),I\left(s\right),V\left(s\right),T\left(s\right)$, for $s\in [-\Delta ,0]$. Therefore, we must guarantee the positivity and boundedness of the solution $\left(I_E\left(t\right),I\left(t\right),V\left(t\right),T\left(t\right)\right)$ of the system (6) at $[0,\infty )$. The following results characterizes these properties.

Theorem 1.

If the initial conditions $I_E\left(0\right)=I_{E_0},I\left(0\right)=I_0,V\left(0\right)=V_0,T\left(0\right)=T_0$ of the mathematical model (6) are positive, then the solutions $\left(I_E\left(t\right),I\left(t\right),V\left(t\right),T\left(t\right)\right)$ of the system (6) are positive for all $t\in [0,\infty ).$

Proof. 

Let us start by noting that the solutions of the differential equations of $I_E\left(t\right)$ and $I\left(t\right)$ given in (6) can be written as

$$\begin{array}{c}{I}_E\left(t\right)=e^{-{\delta }_{I_E}t}\left[I_{E_0}+{\int }_{t-\Delta }^t\beta T\left(s\right)V\left(s\right)e^{{\delta }_{I_E}s}ds\right],\end{array}$$

(8)

$$\begin{array}{c}I\left(t\right)=e^{-{\delta }_{I}t}\left[I_0+{\int }_0^t\beta T(s-\Delta )V(s-\Delta )e^{({\delta }_I-{\delta }_{I_E})s}ds\right].\end{array}$$

(9)

Therefore, the positivity of the solutions $T\left(t\right)$ and $V\left(t\right)$ for all $t>0$, allows us to guarantee the positivity of $I_E\left(t\right)$ and $I\left(t\right)$ and thus of system (6). Thus, for $T\left(t\right)$ given as in (6) we have that $T\left(t\right)>0,$ for all $t\ge 0$. Indeed, suppose that the positivity does not holds, therefore there must be a $t_0>0$ such that $T\left(t_0\right)=0$, $\frac{dT\left(t_0\right)}{dt}\le 0$ and $T\left(t\right)>0$ for all $t\in [0,t_0)$, because the initial condition $T_0>0$. Thus, $T\left(t\right)$ must be negative from some $t_0$. However, in the interval $[0,t_0)$ the function $T\left(t\right)$ is positive, and at point $t_0$ the derivative at $t_0$ is non-positive. Thus, from the fourth equation of model (6), it follows that for $t_0,$

$$\frac{dT\left(t_0\right)}{dt}=\mathsf{\Lambda }-\beta T\left(t_0\right)V\left(t_0\right)-{\mu }_{0}T\left(t_0\right)=\mathsf{\Lambda }>0,$$

which contradicts that $\frac{dT\left(t_0\right)}{dt}\le 0$. Therefore, we must have $T\left(t\right)>0$, for all $t\ge 0$. Now, we affirm that if $V\left(t\right)$ is given by the system (6), it follows that

$$V\left(t\right)>0,\mathrm{for}\mathrm{all}t\ge 0.$$

(10)

Indeed, suppose that there exists a $t_1>0$ such that $V\left(t_1\right)=0$, $\frac{dV\left(t_1\right)}{dt}\le 0$ and $V\left(t\right)>0$ for all $t\in [0,t_1)$. Then it holds that $I\left(t\right)>0$ for all $t\in [0,t_1)$. Otherwise, there should be a $t_2$ such that $0<t_2<t_1$, $I\left(t_2\right)=0$, $\frac{dI\left(t_2\right)}{dt}\le 0$ and $I\left(t\right)>0$ for all $t\in [0,t_2)$. Thus, from the second equation of system (6), if follows that for $t_2\in (\Delta ,t_1)$ it holds

$$\frac{dI\left(t_2\right)}{dt}=\beta T(t_2-\Delta )V(t_2-\Delta )e^{-{\delta }_{I_E}\Delta }-{\delta }_{I}I\left(t_2\right)=\beta T(t_2-\Delta )V(t_2-\Delta )e^{-{\delta }_{I_E}\Delta }.$$

But, $-\Delta <t_2-\Delta <t_1-\Delta <t_1$. From the initial conditions given by (7) and the hypothesis for $V\left(0\right)>0$ it follows that $V(t_2-\Delta )>0$. This, contradicts that $\frac{dI\left(t_2\right)}{dt}\le 0$. Thus, $I\left(t\right)>0$, for all $t\in [0,t_1)$. Next, from third equation of system (6) for $t_1>0,$

$$\frac{dV\left(t_1\right)}{dt}=N{\delta }_{I}I\left(t_1\right)-CV\left(t_1\right)-\beta T\left(t_1\right)V\left(t_1\right)=N{\delta }_{I}I\left(t_1\right)>0.$$

This is a contradiction since $\frac{dV\left(t_1\right)}{dt}\le 0$. Therefore, $V\left(t\right)>0$ for all $t\ge 0$. □

Theorem 2.

The solution $\left(I_E\left(t\right),I\left(t\right),V\left(t\right),T\left(t\right)\right)$ of system (6) is uniformly bounded in $[0,\infty )$.

Proof. 

From system (6), one can see that

$$\frac{d{I}_E\left(t\right)}{dt}+\frac{dI\left(t\right)}{dt}+\frac{dT\left(t\right)}{dt}=\mathsf{\Lambda }-{\delta }_{I_E}{I}_E\left(t\right)-{\delta }_{I}I\left(t\right)-{\mu }_{0}T\left(t\right)\le \mathsf{\Lambda }-M(I_E\left(t\right)+I\left(t\right)+T\left(t\right)),$$

where $M=min\left\{{\delta }_{I_E},{\delta }_I,{\mu }_0\right\}$. This implies that

$$\begin{array}{cc}\hfill I_E\left(t\right)+I\left(t\right)+T\left(t\right)& \le e^{-Mt}\left[I_{E_0}+I_0+T_0+{\int }_0^t\mathsf{\Lambda }{e}^{Ms}ds\right]=e^{-Mt}\left[I_{E_0}+I_0+T_0\right]+\frac{\mathsf{\Lambda }}{M}\left(1-e^{-Mt}\right).\hfill \end{array}$$

Thus,

$$\underset{t\to \infty }{lim\; sup}(I_E\left(t\right)+I\left(t\right)+T\left(t\right))\le \frac{\mathsf{\Lambda }}{M}.$$

Accordingly, $I_E\left(t\right),I\left(t\right)$ and $T\left(t\right)$ are uniformly boundedness. Even more, given $\epsilon >0$, there exists $t_1>0$ such that for all $t\ge t_1$,

$$I\left(t\right)\le \frac{\mathsf{\Lambda }}{M}+\epsilon .$$

Then, since $\frac{dV\left(t\right)}{dt}\le N{\delta }_{I}I\left(t\right)-CV\left(t\right)$, and $V\left(t\right)>0$ for all $t>t_1$, one obtains that

$$\begin{array}{c}\hfill \frac{dV\left(t\right)}{dt}+CV\left(t\right)\le N{\delta }_I\left(\frac{\mathsf{\Lambda }}{M}+\epsilon \right).\end{array}$$

It follows that

$$\begin{array}{c}V\left(t\right)\le V_0{e}^{-Ct}+\left(1-e^{-Ct}\right)\left(\frac{N{\delta }_I}{C}\right)\left(\frac{\mathsf{\Lambda }}{M}+\epsilon \right).\end{array}$$

As $t⟶\infty$, then $V\left(t\right)\le \left(\frac{N{\delta }_I}{C}\right)\left(\frac{\mathsf{\Lambda }}{M}+\epsilon \right).$ Since $ϵ>0$, $V\left(t\right)$ is uniformly boundedness. This completes the proof. □

3.1. Equilibrium Points

The model described by the system of differential Equation (6) has two stationary states, the first one corresponds to the disease-free equilibrium and the second to the endemic equilibrium, which we will denote $P^0$ and $P^{*}$ respectively. To determine both states we must calculate the critical points of the system (6) by setting $\frac{d{I}_E\left(t\right)}{dt}=\frac{dI\left(t\right)}{dt}=\frac{dV\left(t\right)}{dt}=\frac{dT\left(t\right)}{dt}=0$. Thus, we have

$$\begin{array}{cc}\hfill 0& =\beta TV-\beta TV{e}^{-{\delta }_{I_E}\Delta }-{\delta }_{I_E}{I}_E\hfill \\ \hfill 0& =\beta TV{e}^{-{\delta }_{I_E}\Delta }-{\delta }_{I}I\hfill \\ \hfill 0& =N{\delta }_{I}I-CV-\beta TV\hfill \\ \hfill 0& =\mathsf{\Lambda }-\beta TV-{\mu }_{0}T.\hfill \end{array}$$

(11)

The disease-free equilibrium point of a model are solutions of steady state in the absence of infection. For this case, we must consider $I_E=0$, $I=0$, $V=0$, and $T>0$, in the system (11). Then $P^0$ will be of the form $P^0=(0,0,0,T^0)$, where $T^0=\frac{\mathsf{\Lambda }}{{\mu }_0}$. Therefore, $P^0=\left(0,0,0,\frac{\mathsf{\Lambda }}{{\mu }_0}\right).$

On the other hand, we can determine the basic reproductive number using the next generation matrix methodology. With the terms of infection and viral production in the mathematical model (11), matrices F and V are given by

$$\begin{array}{c}\hfill \mathbb{F}=\left(\begin{array}{ccc}0& 0& \beta T^0\left(1-e^{-{\delta }_{I_E}\Delta }\right)\\ 0& 0& \beta T^0{e}^{-{\delta }_{I_E}\Delta }\\ 0& 0& 0\end{array}\right),\mathbb{V}=\left(\begin{array}{ccc}{\delta }_{I_E}& 0& 0\\ 0& {\delta }_I& 0\\ 0& N{\delta }_I& C+\beta T^0\end{array}\right),\end{array}$$

where $T^0=\frac{\mathsf{\Lambda }}{{\mu }_0}$, which it is the number of target cells before infection. Thus, the basic reproductive number ${\mathcal{R}}_0$, is calculated as the spectral radius of the matrix given by

$$\begin{array}{cc}\hfill \mathbb{F}{\mathbb{V}}^{-1}& =\left(\begin{array}{ccc}0& 0& \beta T^0\left(1-e^{-{\delta }_{I_E}\Delta }\right)\\ 0& 0& \beta T^0{e}^{-{\delta }_{I_E}\Delta }\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}\frac{1}{{\delta }_{I_E}}& 0& 0\\ 0& -{\delta }_I& 0\\ 0& -\frac{N}{C+\beta T^0}& \frac{1}{C+\beta T^0}\end{array}\right)\hfill \\ & =\left(\begin{array}{ccc}0& -\frac{N\beta T^0\left(1-e^{-{\delta }_{I_E}\Delta }\right)}{\beta T^0+C}& \frac{\beta T^0\left(1-e^{-{\delta }_{I_E}\Delta }\right)}{\beta T^0+C}\\ 0& -\frac{N\beta T^0{e}^{-{\delta }_{I_E}\Delta }}{\beta T^0+C}& \frac{\beta T^0{e}^{-{\delta }_{I_E}\Delta }}{\beta T^0+C}\\ 0& 0& 0\end{array}\right).\hfill \end{array}$$

Therefore, the basic reproductive number ${\mathcal{R}}_0$ is given by

$$\begin{array}{c}\hfill {\mathcal{R}}_0=\frac{N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{C{\mu }_0+\beta \mathsf{\Lambda }}.\end{array}$$

(12)

Now, the endemic equilibrium point of a model is its steady-state solutions in the presence of infection or disease, for which it must be considered $I_E>0$, $I>0$, $V>0$ and $T>0$ in the system (11.) Then $P^{*}$ will be of the form $P^{*}=(I_E^{*},I^{*},V^{*},T^{*}).$ In this case, from the system (11) the following equalities are obtained

$$I_E^{*}=\frac{\beta T^{*}{V}^{*}(1-e^{-{\delta }_{I_E}\Delta })}{{\delta }_{I_E}}.$$

(13)

$$I^{*}=\frac{\beta T^{*}{V}^{*}{e}^{-{\delta }_{I_E}\Delta }}{{\delta }_I}.$$

(14)

$$I^{*}=\frac{V^{*}}{N{\delta }_I}(C+\beta T^{*}).$$

(15)

$$T^{*}=\frac{\mathsf{\Lambda }}{\beta V^{*}+{\mu }_0}.$$

(16)

Replacing (16) in (14) and (15), is obtains

$$\begin{array}{c}\hfill I^{*}=\frac{\beta \mathsf{\Lambda }{V}^{*}{e}^{-{\delta }_{I_E}\Delta }}{{\delta }_I(\beta V^{*}+{\mu }_0)}=\frac{V^{*}}{N{\delta }_I}\left(C+\frac{\beta \mathsf{\Lambda }}{\beta V^{*}+{\mu }_0}\right).\end{array}$$

(17)

Then

$$\begin{array}{c}\hfill V^{*}=\frac{N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }-C{\mu }_0-\beta \mathsf{\Lambda }}{C\beta }=\frac{\left(C{\mu }_0+\beta \mathsf{\Lambda }\right)\frac{N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{C{\mu }_0+\beta \mathsf{\Lambda }}-1}{C\beta }.\end{array}$$

(18)

Thus, if ${\mathcal{R}}_0>1$, then $N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }-C{\mu }_0-\beta \mathsf{\Lambda }>0$. Hence, we can write $V^{*}$ as

$$\begin{array}{c}\hfill V^{*}=\frac{\left(C{\mu }_0+\beta \mathsf{\Lambda }\right)\left({\mathcal{R}}_0-1\right)}{C\beta }.\end{array}$$

(19)

Next, we replace (18) in (17) to get

$$I^{*}=\frac{N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }-C{\mu }_0-\beta \mathsf{\Lambda }}{{\delta }_I\beta e^{{\delta }_{I_E}\Delta }(N{e}^{-{\delta }_{I_E}\Delta }-1)}=\frac{\left(C{\mu }_0+\beta \mathsf{\Lambda }\right)\left({\mathcal{R}}_0-1\right)}{{\delta }_I\beta e^{{\delta }_{I_E}\Delta }(N{e}^{-{\delta }_{I_E}\Delta }-1)}.$$

(20)

Now, substituting (18) in (16) one gets

$$T^{*}=\frac{\mathsf{\Lambda }C}{\left(C{\mu }_0+\beta \mathsf{\Lambda }\right)\left({\mathcal{R}}_0-1\right)+C{\mu }_0}.$$

(21)

Finally, we replace (18) and (21) in (13) to obtain

$$\begin{array}{cc}\hfill I_E^{*}& =\frac{\left(N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }-C{\mu }_0-\beta \mathsf{\Lambda }\right)\left(e^{{\delta }_{I_E}\Delta }-1\right)}{{\delta }_I\beta e^{{\delta }_{I_E}\Delta }\left(N{e}^{-{\delta }_{I_E}\Delta }-1\right)}=\frac{\left(C{\mu }_0+\beta \mathsf{\Lambda }\right)\left({\mathcal{R}}_0-1\right)\left(e^{{\delta }_{I_E}\Delta }-1\right)}{{\delta }_I\beta e^{{\delta }_{I_E}\Delta }\left(N{e}^{-{\delta }_{I_E}\Delta }-1\right)}\hfill \end{array}$$

(22)

Note that $I_E^{*}>0$, $I^{*}>0$, $V^{*}>0$ and $T^{*}>0$ if only if ${\mathcal{R}}_0=\frac{N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{C{\mu }_0+\beta \mathsf{\Lambda }}>1$. Thus, $N{e}^{-{\delta }_{I_E}\Delta }>\frac{C{\mu }_0+\beta \mathsf{\Lambda }}{\beta \mathsf{\Lambda }}>1$.

3.2. Local Stability of the Equilibrium Points

Theorem 3.

The disease-free equilibrium point $P^0$ of the system (6) is locally asymptotically stable if ${\mathcal{R}}_0<1.$

Proof. 

The eigenvalues of the Jacobian matrix of system (6) evaluated at point $P^0$, are given as the roots of polynomial

$$\begin{array}{c}\hfill (-{\delta }_I-\lambda )(-{\mu }_0-\lambda )\left[{\lambda }^2+\left({\delta }_I+C+\frac{\beta \mathsf{\Lambda }}{{\mu }_0}\right)\lambda +\left({\delta }_{I}C+\frac{{\delta }_I\beta \mathsf{\Lambda }}{{\mu }_0}-\frac{\beta e^{-{\delta }_{I_E}\Delta }N{\delta }_I\mathsf{\Lambda }}{{\mu }_0}\right)\right]=0.\end{array}$$

Therefore, the first two eigenvalues of the Jacobian matrix evaluated at $P^0$ are ${\lambda }_1=-{\delta }_I<0$ and ${\lambda }_2=-{\mu }_0<0.$

Next, since ${\mathcal{R}}_0<1$ then the coefficients of equation

$${\lambda }^2+\left({\delta }_I+C+\frac{\beta \mathsf{\Lambda }}{{\mu }_0}\right)\lambda +\frac{{\delta }_I}{{\mu }_0}\left(C{\mu }_0+\beta \mathsf{\Lambda }\right)\left(1-{\mathcal{R}}_0\right)=0,$$

(23)

are positives. Thus, since there is no sign change between its terms, and by Descartes’ sign rule it is concluded that Equation (23) does not have positive roots. Now, if $\lambda$ is replaced by $-\lambda$ in (23) then

$${\lambda }^2-\left({\delta }_I+C+\frac{\beta \mathsf{\Lambda }}{{\mu }_0}\right)\lambda +\frac{{\delta }_I}{{\mu }_0}\left(C{\mu }_0+\beta \mathsf{\Lambda }\right)\left(1-{\mathcal{R}}_0\right)=0.$$

(24)

Thus, if ${\mathcal{R}}_0<1$ Equation (24) has two sign changes in its terms, and by Descartes’ sign rule it is concluded that there are exactly two negative roots in Equation (23). Therefore, $P^0=\left(0,0,0,\frac{\mathsf{\Lambda }}{{\mu }_0}\right)$ is locally asymptotically stable if ${\mathcal{R}}_0<1.$ □

Theorem 4.

The $P^{*}$ endemic point of the system (6) is locally asymptotically stable if ${\mathcal{R}}_0>1$.

Proof. 

We note that $R_1=N{e}^{-{\delta }_{I_E}\Delta }-1>0$. Thus, the characteristic equation is given by

$$(-{\delta }_{I_E}-\lambda )\left|\begin{array}{ccc}-{\delta }_I-\lambda & \frac{C{e}^{-{\delta }_{I_E}\Delta }}{R_1}& \frac{\beta \mathsf{\Lambda }{R}_1{e}^{-{\delta }_{I_E}\Delta }}{C}-{\mu }_0{e}^{-{\delta }_{I_E}\Delta }\\ \\ N{\delta }_I& -C-\frac{C}{R_1}-\lambda & {\mu }_0-\frac{\beta \mathsf{\Lambda }{R}_1}{C}\\ \\ 0& -\frac{C}{R_1}& -\frac{\beta \mathsf{\Lambda }{R}_1}{C}-\lambda \end{array}\right|=0$$

Therefore, the first eigenvalue of the Jacobian matrix will be ${\lambda }_1=-{\delta }_{I_E}<0$, and the other, are the roots of polynomial

$$\begin{array}{c}\hfill {\lambda }^3+\left({\delta }_I+C+\frac{C}{R_1}+\frac{\beta \mathsf{\Lambda }{R}_1}{C}\right){\lambda }^2+\left(\frac{{\delta }_I\beta \mathsf{\Lambda }{R}_1}{C}+\beta \mathsf{\Lambda }{R}_1+\frac{C{\mu }_0}{R_1}\right)\lambda +{\delta }_I\left(C{\mu }_0+\beta \mathsf{\Lambda }\right)\left({\mathcal{R}}_0-1\right)=0.\end{array}$$

(25)

If ${\mathcal{R}}_0>1$ all the coefficients of Equation (25) are positive, i.e., there is no sign change between their terms, and by Descartes’s sign rule it is concluded that there are no positive roots. Now if $\lambda$ is replaced by $-\lambda$ in (25) it gives us that

$$-{\lambda }^3+\left({\delta }_I+C+\frac{C}{R_1}+\frac{\beta \mathsf{\Lambda }{R}_1}{C}\right){\lambda }^2-\left(\frac{{\delta }_I\beta \mathsf{\Lambda }{R}_1}{C}+\beta \mathsf{\Lambda }{R}_1+\frac{C{\mu }_0}{R_1}\right)\lambda +{\delta }_I\left(C{\mu }_0+\beta \mathsf{\Lambda }\right)\left({\mathcal{R}}_0-1\right)=0$$

(26)

Then, if ${\mathcal{R}}_0>1$ the polynomial (26) has three sign changes between its terms, and by Descartes’s sign rule it is concluded that there are exactly three negative roots of Equation (25). Thus, $P^{*}$ is locally asymptotically stable if ${\mathcal{R}}_0>1.$ □

3.3. Global Stability Analysis of the Mathematical Model

Since the variable $I_E$ does not appear in the other three equations, without loss of generality we will only consider the following three-dimensional system,

$$\begin{array}{cc}\hfill \frac{dI\left(t\right)}{dt}& =\beta T(t-\Delta )V(t-\Delta )e^{-{\delta }_{I_E}\Delta }-{\delta }_{I}I\left(t\right)\hfill \\ \hfill \frac{dV\left(t\right)}{dt}& =N{\delta }_{I}I\left(t\right)-(C+\beta T\left(t\right))V\left(t\right)\hfill \\ \hfill \frac{dT\left(t\right)}{dt}& =\mathsf{\Lambda }-\beta T\left(t\right)V\left(t\right)-{\mu }_{0}T\left(t\right).\hfill \end{array}$$

(27)

To analyze the global stability of the equilibrium points of the system (27), we use the method of the Lyapunov’s functions, and we will use the Volterra function

$$\begin{array}{c}\hfill G\left(x\right)=x-1-lnx\end{array}$$

(28)

for $x>0,$ which is no negative for any $x>0$ and $G\left(x\right)=0$ if and only if $x=1.$

Theorem 5.

If ${\mathcal{R}}_0<1$ then the disease-free equilibrium point $P_f=\left(0,0,\frac{\mathsf{\Lambda }}{{\mu }_0}\right)$ of system (27) is globally asymptotically stable.

Proof. 

We define the Lyapunov functional

$$\begin{array}{c}\hfill \mathcal{V}\left(t\right)=e^{-{\delta }_{I_E}\Delta }{T}^{0}G\left(\frac{T\left(t\right)}{T^0}\right)+\frac{1}{N}V\left(t\right)+I\left(t\right)+e^{-{\delta }_{I_E}\Delta }{\int }_0^{\Delta }\beta T(t-\theta )V(t-\theta )d\theta .\end{array}$$

Now, calculating the time derivative of $\mathcal{V}\left(t\right)$ along the solution of model (27), one gets that

$$\begin{array}{cc}\hfill \frac{d\mathcal{V}\left(t\right)}{dt}=& e^{-{\delta }_{I_E}\Delta }\frac{T\left(t\right)-T^0}{T\left(t\right)}\frac{dT\left(t\right)}{dt}+\frac{1}{N}\frac{dV\left(t\right)}{dt}+\frac{dI\left(t\right)}{dt}+e^{-{\delta }_{I_E}\Delta }\frac{d}{dt}{\int }_0^{\Delta }\beta T(t-\theta )V(t-\theta )d\theta \hfill \\ \hfill =& -{\mu }_0{e}^{-{\delta }_{I_E}\Delta }\frac{{\left(T\left(t\right)-T^0\right)}^2}{T\left(t\right)}-e^{-{\delta }_{I_E}\Delta }\beta T\left(t\right)V\left(t\right)+e^{-{\delta }_{I_E}\Delta }\beta V\left(t\right)T^0+{\delta }_{I}I\left(t\right)-\frac{C}{N}V\left(t\right)-\frac{\beta T\left(t\right)V\left(t\right)}{N}\hfill \\ & +\beta T(t-\Delta )V(t-\Delta )e^{-{\delta }_{I_E}\Delta }-{\delta }_{I}I\left(t\right)-\beta T(t-\Delta )V(t-\Delta )e^{-{\delta }_{I_E}\Delta }-{\delta }_{I}I\left(t\right)+e^{-{\delta }_{I_E}\Delta }\beta T\left(t\right)V\left(t\right)\hfill \\ \hfill \le & -{\mu }_0{e}^{-{\delta }_{I_E}\Delta }\frac{{\left(T\left(t\right)-T^0\right)}^2}{T\left(t\right)}+\frac{C}{N}\left({\mathcal{R}}_0-1\right).\hfill \end{array}$$

Thus, $\frac{d\mathcal{V}\left(t\right)}{dt}<0$ when ${\mathcal{R}}_0<1.$ Therefore, by Lyapunov–LaSalle Invariance Principle, the infection-free equilibrium $E_f$ is globally asymptotically stable. □

4. Design of a $NSFD$ Scheme for the Mathematical Model

The use of differential equations in the modeling of the transmission of infectious diseases has represented a versatile tool to understand better the dynamics of a variety of infectious diseases [7,59,60,74,75,76]. Mathematical models based on differential equations have been useful to study how to reduce the burden of infectious diseases. The models allow the determination of optimal controls and estimate the impact of a variety of virus on the disease dynamics [67,75,77]. One advantage of mathematical models is that different simulations can be performed, and this allows us to analyze different main driven factors of epidemics under a variety of complex scenarios [8,11,59]. However, there are no general formulas that allow the obtaining of precise analytical solutions for many differential equation systems. These solutions exist only occasionally and are often difficult to find, so good approximations are necessary that preserve the qualitative properties of said solution, for which numerical methods have been used, see [24,25,31,48,49,50,78,79,80,81,82,83].

Discrete epidemic models generated by numerical methods contain additional parameters to those that already exist in differential equations, such as the time and space steps. Variations in these additional parameters can generate solutions to the discrete equations that do not correspond to any solution of the original differential equations, producing fictitious bifurcations, artificial chaos, spurious solutions, and false stable states [24,25,26,45,83,84]. Therefore, we must choose numerical discrete schemes that guarantee the qualitative properties of the mathematical models. There are several methods that can be used to obtain accurate solutions. For instance, the Richardson extrapolation on uniform and nonuniform grids or $NSFD$ schemes have been used for that end [45,85,86,87,88,89,90].

Another, important aspect where a robust numerical method plays an important role is when solving inverse problems to estimate the parameters of the model [48,51,56,57,91,92]. Thus, for mathematical models of a variety of virus is of paramount importance to have a robust and efficient numerical method for solving the differential equations [25,49,51]. Usually, when a differential equation is solved numerically a certain tolerance is prescribed and this has an impact in the success in estimating the parameters [48,49,93,94]. In this paper, we deal with a mathematical model that is based on system of differential equations with discrete time delay [17,32,60,72,95,96,97]. There are different numerical methods to deal with this type of equations, and some are analogous to the ones used for ordinary differential equations but with additional issues [50,98,99,100,101,102]. One particular numerical method that we are interested in is by using $NSFD$ schemes to guarantee some properties of the continuous mathematical model. Some previous works using this methodology have been developed for linear and nonlinear delay differential equations [103,104,105,106,107].

For the construction of a discrete numerical scheme that allows us to efficiently approximate the solutions of the system (6), we use the methodology proposed by Ronald Mickens, see [24,25,26,27]. In that order of ideas, for the discrete approximation of the time derivative of a function $X\left(t\right)\in C^1\left(\mathbb{R}\right),$ we define the non-standard derivative as

$$\begin{array}{c}\hfill \frac{d_{N}X\left(t\right)}{dt}=\frac{X(t+h)-X\left(t\right)}{\phi \left(h\right)}+\mathcal{O}\left(h\right),\end{array}$$

(29)

where $\phi \left(h\right)$ is a real positive valued function that satisfies $\phi \left(h\right)=h+\mathcal{O}\left(h^2\right)$, and N is to denote the non-standard derivative.

Although there is no general algorithm for constructing an $NSFD$ schema that approximates the solutions of a given system of differential equations, the following general rules are often useful to correctly design these schemes.

Rule 1.

The discrete derivatives in a numerical scheme must be of the same orders as the continuous derivatives that appear in the differential equation.

Rule 2.

Discrete derivatives may have non-trivial denominators.

Rule 3.

Nonlinear terms that appear in differential equations must have non-local representations.

Rule 4.

The numerical solution must preserve all the special conditions that hold for the solutions of the corresponding differential equations.

Rule 5.

The scheme must not introduce unnecessary or false solutions, i.e., convergence to false steady states.

Let us denote by $I_E^n$, $I^n$, $V^n$ and $T^n$ the approximations of $I_E\left(nh\right)$, $I\left(nh\right)$, $V\left(nh\right)$ and $T\left(nh\right)$, respectively, for $n=0,1,2\dots ,$ and for h size step in time of the scheme. The value of $I_E^{n+1}$ for $n=0,1,\cdots ,$ is calculated using Equation (8) and a quadrature formula. For this case, we use

$$\begin{array}{c}\hfill I_E^{n+1}=e^{-{\delta }_{I_E}\Delta (n+1)h}\left[I_E^n+\frac{m_{1}h}{2}\left(\beta T^{n+1}{V}^{n+1}{e}^{{\delta }_{I_E}\Delta (n+1)h}+\beta T^{n+1-m_1}{V}^{n+1-m_1}{e}^{{\delta }_{I_E}\Delta (n+1-m_1)h}\right)\right],\end{array}$$

(30)

with $\Delta =m_{1}h.$

Next, we make the following non-local approximations of the terms on the right side of the system (27)

$$\left\{\begin{array}{ccc}\hfill \beta T\left(t\right)V\left(t\right)& \to & \beta T^n{V}^n\hfill \\ \hfill -\beta T(t-\Delta )V(t-\Delta )e^{-{\delta }_{I_E}\Delta }& \to & -\beta T^{n-m_1+1}{V}^{n-m_1}{e}^{-{\delta }_{I_E}\Delta }\hfill \\ \hfill -{\delta }_{I_E}{I}_E\left(t\right)& \to & -{\delta }_{I_E}{I}_E^{n+1}\hfill \\ \hfill -{\delta }_{I}I\left(t\right)& \to & -{\delta }_I{I}^{n+1}\hfill \\ \hfill N{\delta }_{I}I\left(t\right)& \to & N{\delta }_I{I}^{n+1}\hfill \\ \hfill -(C+\beta T(t\left)\right)V\left(t\right)& \to & -(C+\beta T^n)V^{n+1}\hfill \\ \hfill -\beta T\left(t\right)V\left(t\right)& \to & -\beta T^{n+1}{V}^n\hfill \\ \hfill -{\mu }_{0}T\left(t\right)& \to & -{\mu }_0{T}^{n+1}\hfill \end{array}\right.$$

(31)

Then, we approximate the derivatives on the left side of the system (27) as follows

$$\left\{\begin{array}{ccc}\hfill \frac{d_{N}I\left(t\right)}{dt}& \to & \frac{I^{n+1}-I^n}{\phi \left(h\right)}\hfill \\ \hfill \frac{d_{N}V\left(t\right)}{dt}& \to & \frac{V^{n+1}-V^n}{\phi \left(h\right)}\hfill \\ \hfill \frac{d_{N}T\left(t\right)}{dt}& \to & \frac{T^{n+1}-T^n}{\phi \left(h\right)}\hfill \end{array}\right.$$

(32)

Consequently, the system (27) can be discretized as an implicit $NSFD$ scheme given by

$$\begin{array}{cc}\hfill \frac{T^{n+1}-T^n}{\phi \left(h\right)}& =\mathsf{\Lambda }-\beta T^{n+1}{V}^n-{\mu }_0{T}^{n+1},\hfill \\ \hfill \frac{I^{n+1}-I^n}{\phi \left(h\right)}& =\beta T^{n-m_1+1}{V}^{n-m_1}{e}^{-{\delta }_{I_E}\Delta }-{\delta }_I{I}^{n+1},\hfill \\ \hfill \frac{V^{n+1}-V^n}{\phi \left(h\right)}& =N{\delta }_I{I}^{n+1}-C{V}^{n+1}-\beta T^n{V}^{n+1}.\hfill \end{array}$$

(33)

And the explicit form is given by

$$\begin{array}{cc}\hfill T^{n+1}& =\frac{\phi \left(h\right)\mathsf{\Lambda }+T^n}{1+\phi \left(h\right)(\beta V^n+{\mu }_0)},\hfill \\ \hfill I^{n+1}& =\frac{\phi \left(h\right)\beta T^{n-m_1+1}{V}^{n-m_1}{e}^{-{\delta }_{I_E}\Delta }+I^n}{1+\phi \left(h\right){\delta }_I},\hfill \\ \hfill V^{n+1}& =\frac{\phi \left(h\right)N{\delta }_I{I}^{n+1}+V^n}{1+\phi \left(h\right)\left(C+\beta T^n\right)},\hfill \end{array}$$

(34)

where $m_1=\frac{\Delta }{h}\in \mathbb{N}$. The initial conditions of scheme (34) are given by

$$\begin{array}{c}\hfill T^j={\xi }_1^j,V^j={\xi }_2^j,j=-m_1,-m_1+1,\cdots ,0.\end{array}$$

The positivity of scheme (34) is trivially satisfied, since for $n>0$ it holds that $T^n,I^n,V^n$ are positive.

Theorem 6.

Let $(T^n,I^n,V^n)$ be a solution of system (34). Then is uniformly bounded for all $n>0.$

Proof. 

From the first equation of scheme (33), one gets that

$$\begin{array}{c}\hfill \frac{T^{n+1}-T^n}{\phi \left(h\right)}=\mathsf{\Lambda }-\beta T^{n+1}{V}^n-{\mu }_0{T}^{n+1}\le \mathsf{\Lambda }-{\mu }_0{T}^{n+1}.\end{array}$$

When $n\to \infty$ and since $\phi \left(h\right)=h+O\left(h^2\right),$ then $\phi \left(h\right)$ coincide with 0 in the limit as $h\to 0$. This implies that $\underset{n\to \infty }{lim\; sup}{T}^n\le \frac{\mathsf{\Lambda }}{{\mu }_0}.$

Next, let ${\mathcal{W}}^n=T^{n-m_1}+e^{{\delta }_{I_E}\Delta }{I}^n.$ From first and second equation of system (33), one obtains

$$\begin{array}{cc}\hfill \frac{{\mathcal{W}}^{n+1}-{\mathcal{W}}^n}{\phi }=& \frac{T^{n-m_1+1}-T^{n-m_1}}{\phi \left(h\right)}+e^{{\delta }_{I_E}\Delta }\frac{I^{n+1}-I^n}{\phi \left(h\right)}=\mathsf{\Lambda }-{\mu }_0{T}^{n-m_1+1}-{\delta }_I{e}^{{\delta }_{I_E}\Delta }{I}^{n+1}\hfill \\ \hfill \le & \mathsf{\Lambda }-d{\mathcal{W}}^{n+1},\hfill \end{array}$$

where $d=\min \left\{{\mu }_0,{\delta }_I\right\}.$ Thus, $\underset{n\to \infty }{lim\; sup}{\mathcal{W}}^n\le \frac{\mathsf{\Lambda }}{d}$. Therefore, $\underset{n\to \infty }{lim\; sup}{I}^n\le \frac{\mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{d}$. Then, give $ϵ>0$ we can choose a $M\in \mathbb{N}$ such that $I^n\le \frac{\mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{d}+ϵ$ for $n\ge M$. Using the last equation of (33) it is concluded that

$$\begin{array}{c}\hfill \frac{V^{n+1}-V^n}{\phi \left(h\right)}\le N{\delta }_I{I}^{n+1}-C{V}^{n+1}\le N{\delta }_I\left(\frac{\mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{d}+ϵ\right)-C{V}^{n+1}\end{array}$$

for $n\ge M+1.$ Then $\underset{n\to \infty }{lim\; sup}{V}^n\le \frac{N{\delta }_I}{C}\left(\frac{\mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{d}+ϵ\right)$, and as is for all $ϵ>0$ it follows that $\underset{n\to \infty }{lim\; sup}{V}^n\le \frac{N{\delta }_I}{C}\left(\frac{\mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{d}\right)$. This completes the proof. Moreover, from Equation (27) it follows that $I_E^n$ is bounded. □

4.1. Equilibrium Points of the $NSFD$ Numerical Scheme

The equilibrium points of the scheme (34) are given by analyzing the behavior of system when n approaches to infinity. Thus, after a few calculations we find that

$$\begin{array}{cc}\hfill I^{*}& =\frac{\phi \left(h\right)\beta T^{*}{V}^{*}{e}^{-{\delta }_{I_E}\Delta }+I^{*}}{1+\phi \left(h\right){\delta }_I}\hfill \\ \hfill V^{*}& =\frac{\phi \left(h\right)N{\delta }_I{I}^{*}+V^{*}}{1+\phi \left(h\right)(C+\beta T^{*})}\hfill \\ \hfill T^{*}& =\frac{\phi \left(h\right)\mathsf{\Lambda }+T^{*}}{1+\phi \left(h\right)(\beta V^{*}+{\mu }_0)}.\hfill \end{array}$$

(35)

Note that the equations of the scheme (35) correspond to Equations (14)–(16). Thus, the critical points of the discrete scheme will coincide in the limit $h\to 0$, with those of the continuous model.

4.2. Local Stability of the $NSFD$ Numerical Scheme

For the study of the local stability of the critical points of the numerical scheme (34) it is necessary to use the following lemma:

Lemma 1.

The roots of the quadratic polynomial ${\lambda }^2-a_1\lambda +a_2=0$, satisfy $\left|{\lambda }_i\right|<1$ for $i=1,2$ if and only if the following conditions hold:

i.

$1-a_1+a_2>0,$

ii.

$1+a_1+a_2>0,$

iii.

$a_2<1.$

Proof. 

See [7]. □

Theorem 7.

The disease-free equilibrium point $P_f=\left(0,0,\frac{\mathsf{\Lambda }}{{\mu }_0}\right)$ of the scheme (34) is locally asymptotically stable if ${\mathcal{R}}_0<1$.

Proof. 

Calculating the eigenvalues of the Jacobian matrix of the system (34) at the disease-free point, we obtain the following characteristic polynomial

$$\begin{array}{c}\hfill \left(\frac{1}{1+\phi \left(h\right){\mu }_0}-\lambda \right)\left|\begin{array}{cc}\frac{1}{1+\phi \left(h\right){\delta }_I}-\lambda & \frac{\phi \left(h\right)\beta e^{-{\delta }_{I_E}\Delta }\mathsf{\Lambda }}{{\mu }_0(1+\phi \left(h\right){\delta }_I)}\\ \frac{{\mu }_{0}hN{\delta }_I}{{\mu }_0+{\mu }_{0}hC+h\beta \mathsf{\Lambda }}& \frac{{\mu }_0}{{\mu }_0+{\mu }_{0}hC+h\beta \mathsf{\Lambda }}-\lambda \end{array}\right|=0.\end{array}$$

Thus, the first eigenvalue is

$$\begin{array}{c}\hfill {\lambda }_1=\frac{1}{1+\phi \left(h\right){\mu }_0}<1.\end{array}$$

The other two eigenvalues, are the roots of quadratic polynomial

$$\begin{array}{c}\hfill {\lambda }^2-\left(\frac{1}{p}+\frac{{\mu }_0}{q}\right)\lambda +\frac{{\mu }_0-\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{{\mu }_{0}pq}=0,\end{array}$$

(36)

where $p=1+\phi \left(h\right){\delta }_I>0$ and $q={\mu }_0+{\mu }_0\phi \left(h\right)C+\phi \left(h\right)\beta \mathsf{\Lambda }$. Next, let $a_1=\frac{1}{p}+\frac{{\mu }_0}{q}$ and $a_2=\frac{{\mu }_0-\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{{\mu }_{0}pq}.$ We have the following affirmations.

1.

If $1>{\mathcal{R}}_0$, it follows that $\phi {\left(h\right)}^2{\delta }_I(C{\mu }_0+\beta \mathsf{\Lambda })>\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }.$ Thus,

$$\begin{array}{c}\hfill {\delta }_I({\mu }_0+{\mu }_{0}hC+h\beta \mathsf{\Lambda })>{\mu }_0{\delta }_I+\phi \left(h\right){\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }⟺\\ \hfill (1+\phi \left(h\right){\delta }_I)q>q+\phi \left(h\right){\mu }_0{\delta }_I+\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }⟺\\ \hfill pq+{\mu }_0>q+\phi \left(h\right){\mu }_0{\delta }_I+\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }+{\mu }_0⟺\\ \hfill pq+{\mu }_0>q+\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }+p{\mu }_0⟺\\ \hfill 1-\left(\frac{1}{p}+\frac{{\mu }_0}{q}\right)+\frac{{\mu }_0-h^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{pq}>0\end{array}$$

Therefore, one gets that $1+a_2>a_1.$

2.

Since $a_1>0$, it is sufficient to prove that $1+a_2>0$. By hypothesis $1>R_0$, then

$$\begin{array}{c}\hfill {\mu }_0+q+\phi \left(h\right){\delta }_I{\mu }_0+\phi {\left(h\right)}^2{\delta }_I\left({\mu }_{0}C+\beta \mathsf{\Lambda }\right)>\phi {\left(h\right)}^2{\delta }_I\left({\mu }_{0}C+\beta \mathsf{\Lambda }\right)>\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }.\end{array}$$

Accordingly

$$\begin{array}{c}\hfill {\mu }_0+q+\phi \left(h\right){\delta }_I({\mu }_0+{\mu }_0\phi \left(h\right)C+\phi \left(h\right)\beta \mathsf{\Lambda })>\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }⟺\\ \hfill {\mu }_0+q+\phi \left(h\right){\delta }_{I}q>\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }⟺\\ \hfill {\mu }_0+(1+\phi \left(h\right){\delta }_I)q>\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }⟺\\ \hfill 1+\frac{{\mu }_0-\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }}{pq}>0\end{array}$$

3.

Given that

$$\begin{array}{c}\hfill {\mu }_0-\phi {\left(h\right)}^2{\delta }_{I}N\beta \mathsf{\Lambda }{e}^{-{\delta }_{I_E}\Delta }<{\mu }_0+{\mu }_0\phi \left(h\right)C+\phi \left(h\right)\beta \mathsf{\Lambda }+\phi \left(h\right){\delta }_{I}q,\end{array}$$

then

$$\begin{array}{cc}\hfill \frac{{\mu }_0-\phi {\left(h\right)}^2\beta e^{-{\delta }_{I_E}\Delta }\mathsf{\Lambda }N{\delta }_I}{pq}& <1,\hfill \end{array}$$

that is $a_2<1$.

Next, by virtue of Lemma 1 we have that the eigenvalues of the Jacobian matrix of the system (35) evaluated at $P_f$ satisfy $\left|{\lambda }_i\right|<1.$ for $i=1,2.$ Then $P_f$ is locally asymptotically stable if ${\mathcal{R}}_0<1$. □

4.3. Global Stability of the $NSFD$ Numerical Scheme

Several authors have used discrete Lyapunov functions to study the global behavior of numerical solutions generated by non-standard finite difference schemes ($NSFD$), see [19,106,107,108,109,110]. For the study of the global stability of the critical points of the numerical scheme (34) it is necessary to use the Lyapunov functions and Equation (28).

Theorem 8.

The disease-free equilibrium point $P_f=\left(0,0,\frac{\mathsf{\Lambda }}{{\mu }_0}\right)$ of the scheme (34) is globally asymptotically stable if ${\mathcal{R}}_0\le 1$.

Proof. 

Let the following Lyapunov function be

$$\begin{array}{c}\hfill {\mathcal{L}}^n=\frac{T^{0}g\left(\frac{T^n}{T^0}\right)+e^{{\delta }_{I_E}\Delta }{I}^n+\beta T^0\left(\frac{1}{C}+\phi \left(h\right)\right)V^n}{\phi \left(h\right)}+\sum _{i=n-m_1}^{n-1}\beta T^{i+1}{V}^i.\end{array}$$

(37)

Using the inequality $\ln z\le z-1,$ the difference of ${\mathcal{L}}^n$ in (37) satisfies

$$\begin{array}{cc}\hfill {\mathcal{L}}^{n+1}-{\mathcal{L}}^n& \le \left(\frac{T^{n+1}-T^0}{T^{n+1}}\right)\left(\mathsf{\Lambda }-\beta T^{n+1}{V}^n-{\mu }_0{T}^{n+1}\right)\hfill \\ & +e^{{\delta }_{I_E}\Delta }\left(\beta T^{n-m_1+1}{V}^{n-m_1}{e}^{-{\delta }_{I_E}\Delta }-{\delta }_I{I}^{n+1}\right)\hfill \\ & +\frac{\beta T^0}{C}\left(N{\delta }_I{I}^{n+1}-C{V}^{n+1}-\beta T^n{V}^{n+1}\right)+\beta T^0\left(V^{n+1}-V^n\right)\hfill \\ & +\beta T^{n+1}{V}^n-\beta T^{n-m_1+1}{V}^{n-m_1}\hfill \\ & \le -{\mu }_0\frac{T^{n+1}-T^0}{T^{n+1}}-e^{{\delta }_{I_E}\Delta }{\delta }_I{I}^{n+1}\left(1-{\mathcal{R}}_0\right).\hfill \end{array}$$

It follows that ${\mathcal{L}}^{n+1}-{\mathcal{L}}^n\le 0$ for all $n\le 0$ if ${\mathcal{R}}_0\le 1$. This means that ${\mathcal{L}}^n$ is monotone decreasing sequence and since ${\mathcal{L}}^n\ge 0$ for $n\ge 0$ then there exists a limit, i.e., $\underset{n\to \infty }{\lim }{\mathcal{L}}^n\ge 0$. Hence $\underset{n\to \infty }{\lim }\left({\mathcal{L}}^{n+1}-{\mathcal{L}}^n\right)=0$, which implies that $\underset{n\to \infty }{\lim }{T}^n=T^0$, $\underset{n\to \infty }{\lim }{I}^n=0$, and from (34) $\underset{n\to \infty }{\lim }{V}^n=0$. By the previous analysis, we conclude that $P_f$ is globally asymptotically stable for scheme (34). □

5. Numerical Simulations Using the NSFD Scheme

In this section, we present some numerical solutions of the mathematical model of HIV. To carry out the simulations we use the constructed $NSFD$ numerical scheme (34). We choose the parameter values based on existing experimental data and previous model studies, see [54,106,111,112]. The values of these parameters are given in

Table 1. Parameter values for the numerical simulations when ${\mathcal{R}}_0<1.$

Parameters	Value
β	4.8 × 10 −6 mm3 day−1
C	3 mm3 day−1
Λ	2.3 day−1
δI	0.24 day−1
${\delta }_{I_E}$	0.05 day−1
μ0	0.0046 day−1
Δ	0.4 day
N	500
V0	1 mm−3
T0	100 mm−3
$I_{E_0}$	0
I0	0

and

Table 2. Parameter values for the numerical simulations when ${\mathcal{R}}_0>1.$

Parameters	Value
β	4.8 × 10 −6 mm3 day−1
C	2.4 mm3 day−1
Λ	23 day−1
δI	0.2 day−1
${\delta }_{I_E}$	0.05 day−1
μ0	0.0046 day−1
Δ	0.4 day
N	500
V0	1 mm−3
T0	100 mm−3
$I_{E_0}$	1
I0	1

. For the first case we choose values such ${\mathcal{R}}_0<1$ and for the second we have the case ${\mathcal{R}}_0>1.$

For the numerical simulations we use a small step size, $h=0.001$. For the discrete derivatives given in system (33), we have many options for the denominator function $\phi$. We have chosen $\phi \left(h\right)=(1-exp(-hp\left)\right)/p$, where $p=\max \left\{A,\Delta ,{\mu }_0,{\delta }_I,{\delta }_{I_E}\right\}$ are parameters of the model and included in the numerical scheme (33). This particular $\phi$ usually provides better numerical results based on previous articles related to $NSFD$ schemes [113,114]. In addition, this option satisfies the asymptotic relation $\phi \left(h\right)=h+O\left(h^2\right)$, and Rule 1. We performed numerical simulations to show that the solutions obtained by the proposed $NSFD$ scheme and the well-known MATLAB routine dde23 agree. A great advantage of the proposed $NSFD$ numerical scheme (34) is that the computation time is smaller. For the first case, the numerical solution using the $NSFD$ scheme is obtained in approximately 0.083992 s, while the dde23 routine spent 2.573003 s. For the second case, the numerical solution using the $NSFD$ scheme is obtained in approximately 0.176461 s, while the routine dde23 spent 11.181812 s. The results obtained are shown in Figure 2 and Figure 3 respectively. It can be seen that the numerical solution of the proposed $NSFD$ numerical scheme (34) and the one obtained by means of the routine dde23 agree.

6. Conclusions

We proposed a mathematical model based on a set of delay differential equations that describe intracellular HIV infection. The model considers three different subpopulations of cells and the HIV virus. The mathematical model is formulated in such a way that takes into account the time between viral entry into a target cell and the production of new virions. Moreover, this time is included using a discrete time delay. We analyzed the local stability of the infection-free and endemic equilibrium states. By using a suitable Lyapunov functional and the LaSalle invariant principle, it is proved that if the basic reproduction ratio is less than unity, the infection-free equilibrium is globally asymptotically stable. In addition, we designed a non-standard difference scheme that preserves some properties of the continuous model. We prove that the constructed NSFD scheme has the same equilibrium points of the continuous model, and the disease-free equilibrium holds the same stability properties. As required by the constraints of the real phenomena, the solutions given by the numerical scheme satisfy positivity and boundedness. The numerical simulations corroborate that the developed NSFD numerical scheme preserves the properties of the continuous model and presents a robust behavior when working with a variety of parameter values.