A finite-difference discretization preserving the structure of solutions of a diffusive model of type-1 human immunodeficiency virus

We investigate a model of spatio-temporal spreading of human immunodeficiency virus HIV-1. The mathematical model considers the presence of various components in a human tissue, including the uninfected CD4+T cells density, the density of infected CD4+T cells, and the density of free HIV infection particles in the blood. These three components are nonnegative and bounded variables. By expressing the original model in an equivalent exponential form, we propose a positive and bounded discrete model to estimate the solutions of the continuous system. We establish conditions under which the nonnegative and bounded features of the initial-boundary data are preserved under the scheme. Moreover, we show rigorously that the method is a consistent scheme for the differential model under study, with first and second orders of consistency in time and space, respectively. The scheme is an unconditionally stable and convergent technique which has first and second orders of convergence in time and space, respectively. An application to the spatio-temporal dynamics of HIV-1 is presented in this manuscript. For the sake of reproducibility, we provide a computer implementation of our method at the end of this work.


Introduction
In this manuscript, we agree that a, b, and T * are real numbers such that a < b and T * > 0.
We fix the spatial domain B = (a, b) and the space-time domain = B × (0, T * ). The notation is used to denote the closure of the set in the usual topology of R 2 , and we use ∂B to represent the boundary of the set B. In this work, we assume that the functions T, U, V : → R are sufficiently smooth. Meanwhile, the constants β, d, k, δ, γ , c, and N represent nonnegative numbers. Also, we define the functions φ T , φ U , φ V : B → R and ψ T , ψ U , ψ V : ∂B × [0, T] → R. Assume additionally that φ W (x) = ψ W (x, 0) holds for each x ∈ ∂B and W ∈ {T, U, V }. Under these conventions and nomenclature, the model of type-1 human immunodeficiency virus (HIV-1) infection of CD4 + T cells with diffusion is described by the onedimensional problem with initial-boundary conditions: This model is a system with diffusion. The functions T(x, t), U(x, t), and V (x, t) represent the normalized densities of the uninfected CD4 + T cells, infected CD4 + T cells, and the free HIV-1 infection particles in the blood, respectively. The physical meanings of the parameters β, d, κ, δ, γ , c, and N are given in Table 1.
In order to express system (1) in an equivalent form, we suppose that T, U, and V are positive solutions of system (1), and let λ ∈ R + be a free constant. Dividing both sides of each equation of the population system by T(x, t) + λ, U(x, t) + λ, and V (x, t) + λ, respectively, and using the chain rule at the left-hand side of each equation, we obtain the following equivalent system: This equivalent form is employed to propose an exponential-type discretization of the continuous problem under investigation. In particular, we provide a Bhattacharya-type discrete scheme to solve the mathematical model (2). The reason to follow this approach obeys the need to preserve some important features of the relevant solutions of this system and to provide an unconditionally stable and explicit numerical solution for our differential model. It is worth pointing out that the mathematical investigation of the human immunodeficiency virus HIV-1 is an interesting avenue of research. In fact, some works investigate mathematical models to estimate HIV-1 virological failure and establish rigorously the role of lymph node drug penetration [1], the global analysis of the dynamics of predictive systems for intermittent HIV-1 treatment [2], mathematical models of cell-wise spread of HIV-1 which include temporal delays [3], models for patterns of the sexual behavior and their relation with the spread of HIV-1 [4], and the long-term dynamics in mathematical models of HIV-1 with temporal delay in various variants of the drug therapy [5]. Some of these models are based on ordinary differential equations, and their analytical study is followed by simulation experiments which assess the validity of the qualitative results. To that end, various numerical methodologies have been designed and analyzed, like some algorithms for simulating the HIV-1 dynamics at a cellular level [6], stem cells therapy of HIV-1 infections [7], fractional optimal control problems on HIV-1 infection of CD4 + T cells using Legendre spectral collocation [8], HIV-1 cure models with fractional derivatives which possess a nonsingular kernel [9], stochastic HIV-1/AIDS epidemic models in two-sex populations [10], among other reports [5,[11][12][13].
Notice that system (2) is an integer-order diffusive extension of some HIV-1 propagation models available in the literature [9,14]. The use of such a system is due to the current information available of the mechanisms of CD4 + T cells and free HIV-1 infection particles in the blood. In our investigation, we propose a two-level finite difference discretization of (2). Our approach hinges on an exponential-type discretization of the mathematical model, and we prove that the numerical model has various numerical and analytical properties which make it a useful research tool in the study of the propagation of HIV-1. For instance, we prove that the scheme is capable of preserving the positivity and boundedness of solutions. This feature is of the utmost importance in view that the variables under investigation are densities [15]. The properties of consistency, stability, and convergence are thoroughly established in this work. In particular, we show that the scheme is unconditionally stable, and that it has first order of convergence in time and second order in space. We provide some simulations to assess the validity of the theoretical results. Moreover, the computer implementation of the scheme used to obtain the simulations is provided in the Appendix at the end of this work.
Before we begin our study, we must mention that the system under investigation (1) has attracted the attention of these authors due to many important reasons. As we pointed out before, the mathematical model is motivated by various particular models available in the literature which do not consider the presence of diffusion. Those systems are described by ordinary differential equations, whence the investigation of their diffusive generalizations is an important topic of research. Indeed, the consideration of a nonconstant diffusion gives rise to a more realistic and complex scenario. Physically and mathematically, the study of (1) would yield more interesting results. From the numerical perspective, it becomes necessary to possess a reliable tool to investigate the solutions of the mathematical model. After the theoretical and computational investigation of the scheme, researchers in the area would possess a means to obtain trustworthy results to propose predictions on the propagation of HIV-1 in the human body.

Numerical model
In this stage, we introduce discrete operators to provide a discrete model to approximate the analytical solutions of continuous problem (2). Our approach employs finite differences, and the method to solve the continuous problem is introduced herein. The main structural features of the proposed scheme is rigorously established in the second half of this section.
Let us define the sets I q = 1, 2, . . . , q and I q = {0} ∪ I q for each q ∈ N. Let M and K be natural numbers. We define the set ∂J = I M ∩ ∂B and consider discrete partitions corresponding to the intervals [a, b] and [0, T * ] as respectively. In the first partition, the value Let W be any of T, U, or V . We introduce the following discrete quantities for each m ∈ I M-1 and each k ∈ I K-1 : It is well known that the first operator yields a first-order estimate for the partial derivative of W with respect to t at the point (x m , t k ), while the second operator yields a secondorder estimate of the second partial derivative of W with respect to x at the point (x m , t k ).
Substituting these discrete operators at the time t k into model (2), we reach the next finitedifference scheme to estimate the solutions of (2) at (m, k) ∈ I M-1 ×Ī K-1 : such that It is clear that this numerical model is a two-step exponential discretization of the continuous problem (2). Indeed, using the discrete operators, it is an easy algebraic task to check that (7) can be equivalently rewritten as follows: where a k W ,m = h -2 W k m+1 and e k W ,m = h -2 W k m-1 for each m ∈ I M-1 , each k ∈Ī K-1 , and W ∈ {T, U, V }. Notice that each of the three equations in (8) can be expressed as T k+1 Here, the expressions of the functions F T , F U , and F V are as follows: In turn, each function g T , g U , and g V is given by For convenience, we define the (M + 1)-dimensional real vectors for each k ∈Ī K . In general, we say that a vector W ∈ R is positive if all the components are positive. In such a case, we use the notation W > 0. We say that W is bounded from above by s ∈ R if all the components of W are less that s, in which case we employ the notation W < s. Finally, if s is a positive number, then we use 0 < W < s to represent that W > 0 and W < s.
The following results show the existence and uniqueness of the solutions of (8) along with the preservation of the constant solutions.

Theorem 1 (Existence and uniqueness) Let k ∈Ī
Proof The numbers T k m + λ, U k m + λ, and V k m + λ are greater than zero. As a consequence, the real numbers T k+1 are defined uniquely, whence the existence and uniqueness readily follow.

Theorem 2 (Constant solutions)
For each k ∈ I K , let T k , U k , and V k be the zero vectors of dimension M + 1. Then the sequences Proof By the hypothesis, the vectors T 0 = U 0 = V 0 = 0 satisfy the initial-boundary con- The conclusion follows using induction.
The next lemma is an important tool to show the positivity and boundedness of the solutions of the discrete system. The proposition is a result from real analysis, and its proof is established through the mean value theorem.
and some λ ∈ R. Suppose that g and ϕ are differential and that, for each w ∈ [0, 1], the inequality holds. Then F is an increasing function in [0, 1].

Lemma 4
Let λ > 0 and k ∈Ī K-1 . Define the following positive constants: and assume that 0 < T k < 1, 0 < U k < 1, and After some algebra, it is possible to see that where As a consequence, observe that Let W ∈ {T, U, V }. In the following, R φ W and R ψ W represent the ranges of the functions φ W and ψ W , respectively, over the interval [0, 1].

Theorem 5 (Positivity and boundedness)
Let λ > 0, and suppose that the following inequalities are satisfied: Proof We use induction to reach the conclusion. By hypothesis, the conclusion of this theorem is satisfied for k = 0, so let us assume that it holds also for some k ∈Ī K-1 . Lemma 3 assures that the functions F T , F U , and F V are increasing on [0, 1]. Let m ∈ I K-1 . If β = 0 and T k m = U k m = V k m = 0, then it follows that On the other hand, the hypothesis establishes that This and (9) show that F T (1) < 1, F U (1) < 1, and F V (1) < 1. Notice that the functions all are numbers in the set (0, 1) for each m ∈ I M-1 . Using the data at the boundary, we reach 0 < T k+1 m < 1, 0 < U k+1 m < 1, and 0 < V k+1 m < 1. The conclusion follows using induction.
As a conclusion of this section, the numerical methodology is a structure-preserving scheme to approximate the solutions of (2). In this manuscript, the concept of 'structure preservation' or 'dynamical consistency' refers not only to the capacity of discrete models to keep discrete versions of some physical features. In this context, these notions refer also to the capability of a numerical method to be able to conserve some mathematical characteristics of the solutions of interest of continuous paradigms, like positivity [16], boundedness [17,18], monotonicity [19], and convexity of approximations [20], among other physically relevant features [21].

Numerical properties
In this stage, we present the main numerical features of the finite-differences scheme (8). More precisely, we are interested in proving consistency, unconditional stability, and convergence. To show the consistency of the numerical scheme, we require the following continuous operators: for each (x, t) ∈ . Also, we define the difference operators for each m ∈ I M-1 and k ∈Ī K-1 . For the remainder, the symbols · 2 and · ∞ are used to denote the Euclidean and the maximum norms in R M+1 , respectively.

Theorem 6 (Consistency) If T, U, V ∈ C 4,3
x,t ( ) and λ > 0, then there exist positive constants C T , C U , and C V , which are independent of τ and h, such that for each m ∈ I M-1 , k ∈Ī K-1 , and W ∈ {T, U, V }.
Proof We prove the consistency only for W = T, the consistencies for W = U and W = V are proved in a similar fashion. To that end, we employ the usual arguments based on Taylor polynomials. As a consequence of the hypotheses on the regularity of T, there exist positive constants C T,1 and C T,2 independent of τ and h such that, for each (m, k) ∈ I M-1 × I K-1 , On the other hand, observe that Finally, the conclusion follows now from the triangle inequality after we define the positive constant C T = max{C T,1 , C T,2 }, which is independent of τ and h. Similarly, we may show the inequalities corresponding to W = U and W = V .
Under the assumptions of this theorem, there is a positive constant C which is independent of τ and h with the property that, for each m ∈ I M-1 , k ∈Ī K-1 , and W ∈ {T, U, V }, Indeed, observe that C = C T ∨ C U ∨ C V is a constant which satisfies this inequality. This fact will be employed when we prove the convergence of scheme (7). In our next step, we show the stability of the proposed scheme. To that end, we fix two systems of initial-boundary data which are labeled . The corresponding solutions of (7) are represented respectively by (T, U, V ) and ( T, U, V ). In particular, notice that the following result proves that the scheme is unconditionally stable. (7), and let λ > 0. Suppose that the hypotheses of Theorem 5 are satisfied for both triplets (T, U, V ) and ( T, U, V ). There is a constant C, which is independent of the initial data such that

Theorem 7 (Stability) Let and be two sets of initial-boundary conditions for problem
Here, the function (T) is defined by It is readily checked that the function k m is of class C 1 ([0, 1] (M+1) ) for each k ∈ I K . As a consequence of this, the number As consequence, note that, for each m ∈ I M-1 , where Using (51), it is clear that Finally, recursion shows now that the inequality T k - Similarly, we can prove that there exist positive constants C U and C W , with the property Finally, we study the convergence of the numerical scheme (8). In the next result, we let (T, U, V ) be a solution of differential problem (2) associated with the set of initial- Meanwhile, the numerical solution obtained through the discrete model (8) is denoted by ( T, U, V ).

Theorem 8 (Convergence) Let be a set of initial-boundary data which are bounded in
(0, 1), and let λ > 0. Assume that problem (2) has a unique solution bounded in (0, 1) such that T, U, V ∈ C 4,3 x,t ( ). Suppose that the conditions of Theorem 5 hold, and let For each W = T, U, V , there is a constant C W independent of τ and h such that, for each Proof Beforehand, notice that Theorem 5 assures that positive and bounded solutions for discrete problem (7) By Theorem (6), there exists C 0 > 0 such that |R k m | ≤ C 0 (τ + h 2 ) for each m ∈Ī M-1 and k ∈Ī K-1 . Using the definitions of the discrete operators in equations (56) and (57), we have that for each m ∈Ī M-1 and k ∈Ī K-1 . Subtracting T k m of T k m , we have where e k = e k 0 , e k 1 , e k 2 , . . . , for each m ∈Ī M-1 and k ∈Ī K-1 . In addition, we let The constants k m and C k m are as in the proof of Theorem 7. Moreover, all the constants C k m are in the interval [0, 1], therefore C is an element of the interval [0, 1]. Theorem 6 implies now that, for each k ∈Ī K-1 , Taking the sum on both ends of the previous inequality and using the initial data, we have where l ∈Ī K-1 and C T = C 0 DT * . The conclusion of this result has been reached now when W = T. Analogously, we may easily prove the inequality of the conclusion when W = U and W = V .

Application
In this section, we show some computer simulations obtained using the finite-difference scheme (8). Beforehand, notice that the discrete model is an explicit scheme. To describe its computational implementation, for each k ∈Ī K+1 , we redefine the real vectors T k , U k , and V k as follows: These vectors belong to the set R M+1 + , where R + is the system of positive numbers. Also, we define the vectors of initial conditions φ T , φ U , φ V ∈ R M-1 + as follows: Meanwhile, for each k ∈Ī K , we define the vectors ψ k T , ψ k U , ψ k V ∈ R M+1 + of the boundary conditions through With the previous definitions and for each k ∈Ī K , we express the discrete model (8) in vector form as follows: For each k ∈Ī K and W ∈ {T, U, V }, the vectors a k W and e k W are defined as follows: Using the boundary conditions, we readily have that W k M = ψ W (x M , t k ) and W k 0 = ψ W (x 0 , t k ) for each k ∈Ī K+1 and W ∈ {T, U, V }. So, for all k ∈Ī K , the vector form of the finitedifference scheme (8) is defined as follows: The next experiments employ variations of the computer code in the Appendix, which is a basic computational implementation of our finite-difference scheme. The parameter values and the type of initial conditions are motivated by data used in the literature for similar models but without diffusion [9,14].

Figure 1
Snapshots of the approximate solutions T , U, and V of (1) for each row at the time t = 0, t = 1, t = 3, and t = 5. The parameters employed are λ = 0.1, β = 2, d = 5.05, κ = 0.03, δ = 0.016, γ = 3.0, c = 20.4, and N = 1000. We let φ T (x), φ U (x), and φ V (x) be as in Example 9. The approximations were calculated using our implementation of (76) in the Appendix, with τ = 0.00001 and h = 1 Example 9 In this example, we let B = (0, 100) and define the initial conditions φ T , φ U , and φ V as φ T (x) = x 2 , φ U (x) = 0, and φ V (x) = 1 for each x ∈ B. These initial data describe an initial state in which no infected CD4 + T cells are present, and the entire medium is formed from free HIV-1 infection particles. In turn, the normal density of the uninfected CD4 + T cells increases in the linear medium considered herein. Let λ = 0.1, β = 2, d = 5.05, κ = 0.03, δ = 0.016, γ = 3.0, c = 20.4, and N = 1000. Under this situation, Fig. 1 provides snapshots of the normalized solutions T, U, and V at the times t = 0, t = 1, t = 3, and t = 5. The graphs show that the solutions are positive and bounded in accordance with the results obtained in the previous sections. In our simulations, we used the following computational parameters τ = 0.00001 and h = 1.
Example 10 Let B = (0, 100) be as in the previous example, and use the initial conditions φ T (x) = -x 2 + 100x, φ U (x) = 0, and φ V (x) = 1 for each x ∈ B. Set the model and computational parameter values as before. Under these circumstances, Fig. 2  Before closing this section, we would like to point out the biological meaningfulness of the figures obtained in the previous examples. To start with, graphs (a), (d), (g), and (j) of Fig. 1 represent the evolution with respect to the time of the quantity of the uninfected CD4 + T cells. Biologically, the quantity of these cells tends to decrease since there is a significant interaction of the HIV-1 infection with the CD4 + T cells, resulting in an increment of them. The respective increments of the infected CD4 + T cells with respect to the time are shown in graphs (b), (e), (h), and (k). Moreover, the quantity of the infected CD4 + T cells could decrease since an infected cell could die or return to being an uninfected cell. Obviously, this phenomenon is biologically possible. In turn, graphs (c), (f ), (i), and (l) represent the evolution of the HIV-1 infection. From these graphs, it is easy to see that the infection is decreasing with respect to the time due to the presence of an active death rate. The interpretation of the graphs in Fig. 2 is analogous.

Conclusions
In this manuscript, we numerically studied a coupled model consisting of three diffusive nonlinear partial differential equations. The system under investigation is a biological model which describes the interaction of the HIV-1 infection with the CD4 + T cells. One of the equations of the model describes the rate of change of the density of the uninfected CD4 + T cells, the second describes the rate of change of the infected CD4 + T cells, and the third governs the rate of change of the free HIV-1 infection. The differential system was discretized using finite differences to estimate the analytical solutions. The technique that we used in this work is an exponential type that maintains the most important characteristics of the solutions of the continuous model. More concretely, the method was motivated by the well-known family of Bhattacharya exponential-type schemes [22][23][24]. Bhattacharya's discretizations have been employed to derive computational techniques to solve various nonlinear partial differential equations [25][26][27][28]. As it is well known, the main advantage of this family of models is its computational efficiency.
The scheme presented in this work was analyzed to study its most important properties. The most important structural features proved in this work were related to the unique solvability of the discrete model. We also established that the scheme is able to preserve the nonnegativity and the boundedness of the estimations. These properties are highly relevant in light that the functions under investigation represent densities which are positive and bounded. From the numerical point of view, we proved the consistency of the scheme. Moreover, the method is stable, and it converges to the exact solutions with first order in the temporal variable and second order in the spatial. Finally, we provided some computational simulations to illustrate the capability of the scheme to preserve the positivity and the boundedness of the numerical solutions. A Matlab implementation of the method is provided in the Appendix for reproducibility purposes. It is worth pointing out that a study of the mathematical model (1) in two dimensions can be easily performed by extending the theoretical results of this work. Also, an implementation of our scheme in two dimensions is also easily feasible.
Before we close this work, there are various important comments that require to be thoroughly addressed. To start with, it is important to point out that there are some works in the literature where fractional-type models like (1) without diffusion have been investigated [9,14]. Naturally, one would wonder which are the effects of considering a fractional diffusion in such HIV-1 systems. At this point, it is important to mention that one of the authors of the present manuscript has devoted part of his efforts to develop numerical methods for Riesz space-fractional partial differential equations [29,30]. In that context, the differentiation order of the diffusion terms affect the speed of propagation of the spread of effects into the medium. Of course, it would be interesting to propose and analyze numerical models for fractional forms of the system under current investigation. However, the meaningfulness of the use of fractional derivatives in the realistic investigation of HIV-1 may be still questionable. Indeed, not many medical journals employ fractional operators to model the propagation of HIV-1, though the problem is mathematically interesting and challenging.