Comparative Study of Some Numerical Schemes for a Fractional Order Model of HIV Infection Treatment

ABSTRACT

A fractional order mathematical model that already exists in the literature, was considered. This model was established to study the effects of medicinal treatment in people infected with the human immunodeficiency virus (HIV). The importance of this study is that the model evaluates, among other parameters, the density of healthy and HIV-infected CD4⁺ T cells. These data are very necessary for the subject infected by the virus given the effects that an antiretroviral treatment causes in it. The objective of this work is to consider several numerical schemes that involve fractional derivatives in order to compare their behaviors and to obtain a good approximation of the mentioned model solution. Convergence of these schemes will be studied as well as sensitivity with respect to the variation of the parameters η (drug efficacy) and α (fractional derivative order). Furthermore, through the collection of medical records of people living with HIV, it is intended to determine the optimal fractional derivative order for the model and to compare it with the classical model.

Keywords:
Fractional calculus; HIV; numerical schemes

1 INTRODUCTION

In Ferrari et al. ¹³, a fractional order model for the treatment of HIV infection presented in Arafa et al. ³ was considered and properties of the model solution were found. Since this model consists of a coupled system of fractional differential equations, it is not possible to determine the exact solution. Therefore, in this work, several numerical schemes that allow obtaining approximate solutions of it are considered. The behavior of each of them will be studied on this particular model, in order to be able to choose some of the schemes for subsequent analysis. It will be analyzed what happens to the approximate solution when the fractional derivative order approaches to one, comparing with the classical case. We will also see how the coefficient η, which is the efficacy of the reverse transcriptase inhibitor, changes the numerical solution. Once the most appropriate numerical scheme has been chosen, it will be used to study the fractional derivative order that best adapts the model to medical records collected from people living with the virus.

In Section 2, some definitions of the Fractional Calculus used to establish the model of the problem to be studied are given.

In Section 3, the mathematical description of the problem set for the study of the reaction of the immune system of a patient infected with the HIV to the action of the medication received is given.

In Section 4, we describe the numerical schemes to be used to obtain convergence results with respect to the time step. In addition, the behavior of the numerical response of the problem with respect to changes in both the fractional derivative order α and the drug efficacy coefficient η is investigated.

In Section 5, we compare numerical solutions with real data corresponding to patients, who were selected taking into account that a large amount of data were available over time from the start of treatment and that the person had good adherence to treatment.

Finally, in Section 6 we give the corresponding conclusions after analyzing all the information obtained from the computational work.

2 FRACTIONAL CALCULUS

Fractional derivatives have been widely applied in many fields that have seen remarkable growth in the last three decades. For example: heat transfer models that support history, viscoelasticity, electrical circuits, electronic chemistry, economics, physical polymers, and even biology is connected with fractional derivatives. ^{1,16,17,20,21}

As regards HIV, fractional order differential equations are naturally related to memory systems that exist in most biological systems ². Fractional derivatives involve essential characteristics of cellular rheological behavior and have been very successful in the field of rheology ⁷.

For the concept of fractional derivative in the model, Caputo’s definition is adopted, since it has the advantage of properly dealing with initial value problems ⁶.

For a function $f:{\mathrm{ℝ}}^{+}\to \mathrm{ℝ}$, the Riemann-Liouville fractional integral of order α > 0 is given by the expression $J^{\alpha }f\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{\left(t-q\right)}^{\alpha -1}f\left(q\right)dq$ where

$J^{0}f\left(t\right)=f\left(t\right),\alpha >0,t>0.$

The Caputo fractional derivative of order α of a continuous function $f:{\mathrm{ℝ}}^{+}\to \mathrm{ℝ}$ is given by $D_C^{\alpha }f\left(t\right)=J^{m-\alpha }\left(D^{m}f\left(t\right)\right)$, where $m-1<\alpha \le m,m\in \mathrm{\mathbb{N}}$.

3 MATHEMATICAL DESCRIPTION OF THE MODEL

The mentioned model is the one presented in Arafa et al. ³ in its fractional version and in Srivastava et al. ²³ in its classical version. In those, the presence of a reverse transcriptase inhibitor is considered, as it is the most recurrent in treatments.

The variables used in the model are: T (density of susceptible CD4⁺ T cells), I (density of CD4⁺ T cells infected before reverse transcription and therefore do not produce virus), V (density of CD4⁺ T cells infected in those that reverse transcription was completed and are capable of producing viruses), L (virus density).

Then the mathematical problem to be solved is posed as: Find the temporary functions (T, I,V, L), for t ≥ 0, such that they satisfy the following initial value problem:

$\left\{\begin{array}{l}{D}_C^{\alpha }\left(T\right)=s-kLT-\upsilon T+\left(\eta \epsilon +b\right)I\\ D_C^{\alpha }\left(I\right)=kLT-\left({\mu }_1+\epsilon +b\right)I\\ D_C^{\alpha }\left(V\right)=\left(1-\eta \right)\epsilon I-\delta V\\ D_C^{\alpha }\left(L\right)=N\delta V-cL\end{array}\right.$ (3.1)

where 0 < α < 1 is the fractional derivative order, and s, k, µ, η, ε, b, µ1, δ , N, c are positive parameters. The initial conditions that we consider for the problem are $\left(T\left(0\right),I\left(0\right),V\left(0\right),L\left(0\right)\right)=\left(300,10,10,10\right)$. If α = 1 then the classic formulation of the problem is obtained.

In Ferrari et al. ¹³ it is shown that there is one and only one solution $\left(T\left(t\right),I\left(t\right),V\left(t\right),L\left(t\right)\right)$ of the problem modeled in (3.1), where ${\mathrm{ℝ}}_{+}^4$ is a positively invariant domain and that there are two equilibrium points of the problem described in (3.1): a non-infectious status $E_1=\left(\frac{s}{\mu },0,0,0\right)$ and an infectious status $E_2=\left(\overline{T},\overline{I},\overline{V},\overline{L}\right)$, where $\overline{T}=\frac{\left({\mu }_1+\epsilon +b\right)c}{NK\epsilon \left(1-\eta \right)},\overline{I}=\frac{s-\overline{T}\mu }{\epsilon \left(1-\eta \right)+{\mu }_1},\overline{V}=\frac{\left(1-\eta \right)\epsilon \overline{I}}{\delta }$ and $\overline{L}=\frac{N\overline{V}\delta }{c}$.

4 NUMERICAL SCHEMES CONSIDERED

The Adomian Decomposition Method (ADM) and the Variational Iterative Method (VIM) are relatively new approaches that provide an approximate analytical solution to linear and nonlinear problems, and are particularly valuable as tools for scientists and applied mathematicians, since they provide visible and immediate terms of the analytical solutions, as well as approximate numerical solutions of both linear and nonlinear differential equations. Another method that has been developed in recent years is the Homotopy Perturbation Method (HPM). However, all these methods are effective for small times, that is t < 1, and the model studied, which has the t measured in days, requires to know an approximation of the solution for large value of t, for example t = 300 days.

There are also a few other numerical methods for fractional differential equations that have been presented in the mathematical literature ^4,18, however many of them are used for very specific types of differential equations, sometimes only linear equations or even smaller classes. That is why we will start with the simple but no less efficient Generalized Euler and Generalized Trapezoidal methods, which will be developed below.

Considering the function

$$\begin{gathered} \overline{f}\left(t,\left(T\left(t\right),I\left(t\right),V\left(t\right),L\left(t\right)\right)\right) \\ =\left(s-kLT-\mu T+\left(\eta \epsilon +b\right)I,kLT-\left({\mu }_1+\epsilon +b\right)I,\left(I-\eta \right)\epsilon I-\delta V,N\delta V-cL\right), \end{gathered}$$

the problem (3.1) can be rewritten in vector form as:

$\left\{\begin{array}{l}{D}_C^{\alpha }\left(\overline{Y}\right)=\overline{f}\left(t,\overline{Y}\right),\\ \overline{Y}\left(0\right)=Y_1\end{array}\right.$ (4.1)

with $\overline{Y}\left(t\right)=\left(T\left(t\right),I\left(t\right),V\left(t\right),L\left(t\right)\right),t<t_{final}$ the maximum observation time of the patient and $Y_1=\left(T\left(0\right),I\left(0\right),V\left(0\right),L\left(0\right)\right)$ the vector of initial conditions.

In order to obtain a numerical approximation to the solution $\overline{\mathbf{Y}}$ of the initial value problem defined in (4.1), we consider n the number of time steps and we define:

1. $\Delta t=\frac{t_{final}}{n},t_1=0,t_{j+1}=j\Delta t,j=1,\dots ,n-1,$

2. ${\overline{\mathit{Y}}}_{j+1}\approx \overline{Y}\left(t_{j+1}\right),j=1,\dots ,n-1.$

In the first instance, the following methods are used: Generalized Explicit Euler Method (GEEM), Generalized Implicit Euler Method (GIEM) and Generalized Trapezoidal Method (GTM) ^19,22. Given values of α in the interval (0, 1) are considered for these methods.

As it will be seen in the following sections, convergence results were obtained with respect to the time step but in addition we studied how the fractional derivative order α and the η parameter change the solution.

• Generalized Explicit Euler Method (GEEM), which is derived using the Generalized Taylor formula to expand y(t) around t ₀=0. For j = 1,…, n it is defined:

$t_j=j\Delta t,{\overline{\mathbf{Y}}}_j={\overline{\mathbf{Y}}}_{j-1}+\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +1\right)}f\left(t_{j-1},{\overline{\mathbf{Y}}}_{j-1}\right)$ (4.2)

• Generalized Implicit Euler Method (GIEM). For t _j = jΔt, with j = 1,..., n it is defined:

$f\left(t_j,{\overline{\mathbf{Y}}}_j\right)=\frac{1}{\Gamma \left(2-\alpha \right){\left(\Delta t\right)}^{\alpha }}\sum _{k=1}^j\left(k^{1-\alpha }-{\left(k-1\right)}^{1-\alpha }\right)\left[{\overline{\mathbf{Y}}}_{j-k+1}-{\overline{\mathbf{Y}}}_{j-k}\right]$ (4.3)

• Generalized Trapezoidal Method (GTM), which is derived using at each time step the GEEM as predictor and the Generalized Trapezoidal Rule as corrector. For j = 1,..., n is defined j = jΔt and:

$$\begin{gathered} {\overline{\mathbf{Y}}}_{\mathbf{j}}=\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +2\right)}\left({\left(j-1\right)}^{\alpha +1}-\left(j-\alpha -1\right)j^{\alpha }\right)f\left(t_0,{\overline{\mathbf{Y}}}_0\right)+y\left(0\right) \\ +\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +2\right)}\sum _{i=1}^{j-1}\left({\left(j-i+1\right)}^{\alpha +1}-2{\left(j-i\right)}^{\alpha +1}+{\left(j-i-1\right)}^{\alpha +1}\right)f\left(t_i,{\overline{\mathbf{Y}}}_i\right) \\ +\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +2\right)}f\left(t_j,{\overline{\mathbf{Y}}}_{j-1}+\frac{{\left(\Delta t\right)}^{\alpha }}{\Gamma \left(\alpha +1\right)}f\left(t_{j-1},{\overline{\mathbf{Y}}}_{j-1}\right)\right) \end{gathered}$$ (4.4)

In the following, we consider fractional-order numerical schemes under the guidelines described in Changpin et al. ⁵.

• Explicit Fractional Euler Method (EFEM): the approximation of the fractional derivative ${\left[D_{0,t}^{-\alpha }\right]}_{t=t_{n+1}}$ is done by a left-hand rectangular formula:

${\overline{Y}}_{n+1}={\overline{Y}}_1\frac{\Delta t^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{j=1}^n\left[{\left(n-j+1\right)}^{\alpha }-{\left(n-j\right)}^{\alpha }\right]\overline{f}\left(t_j,{\overline{Y}}_j\right),n=1,\dots ,N-1.$ (4.5)

• Implicit Fractional Euler Method (IFEM): in this case, the approximation of the fractional derivative ${\left[D_{0,t}^{-\alpha }\right]}_{t=t_{n+1}}$ is done by a right-hand rectangular formula:

${\overline{Y}}_{n+1}={\overline{Y}}_1+\frac{\Delta t^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{j=1}^n\left[{\left(n-j+1\right)}^{\alpha }-{\left(n-j\right)}^{\alpha }\right]\overline{f}\left(t_{j+1},{\overline{Y}}_{j+1}\right),$ (4.6)

• with n = 1 ,..., N − 1.

• Fractional Weighted Derivative Method (FWDM): for θ ∈ (0, 1) it is proposed:

$$\begin{gathered} {\overline{Y}}_{n+1}={\overline{Y}}_1 \\ +\frac{\Delta t^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{j=1}^n\left[{\left(n-j+1\right)}^{\alpha }-{\left(n-j\right)}^{\alpha }\right]\left[\theta \overline{f}\left(t_j,{\overline{Y}}_j\right)+\left(1-\theta \right)\overline{f}\left(t_{j+1},{\overline{Y}}_{j+1}\right)\right], \end{gathered}$$ (4.7)

• where n = 1 ,..., N − 1.

• Adams Method (AM): it proposes to approximate the fractional derivative ${\left[D_{0,t}^{\alpha }\right]}_{t=t_{n+1}}$ from a trapezoidal formula but, as it is an implicit scheme, it incorporates a predictor-corrector method resulting in the following expression:

${\overline{Y}}_{n+1}={\overline{Y}}_1+\sum _{j=1}^n{a}_{j,n+1}\overline{f}\left(t_j,{\overline{Y}}_j\right)+a_{n+1,n+1}\overline{f}\left(t_{n+1},{\overline{Y}}_{n+1}^P\right),n=1,\dots ,N-1,$ (4.8)

• where

${\overline{Y}}_{n+1}^P={\overline{Y}}_1+\frac{\Delta t^{\alpha }}{\Gamma \left(\alpha +1\right)}\sum _{j=1}^n\left[{\left(n-j+1\right)}^{\alpha }-{\left(n-j\right)}^{\alpha }\right]\overline{f}\left(t_j,{\overline{Y}}_j\right),n=1,\dots ,N-1,$ (4.9)

• and a _{j,n + 1} =

$\frac{\Delta t^{\alpha }}{\Gamma \left(\alpha +2\right)}\left\{\begin{array}{l}{n}^{\alpha +1}-\left(n-\alpha \right){\left(n+\right)}^{\alpha }j=1,\\ {\left(n-j+2\right)}^{\alpha +1}-2{\left(n-j+1\right)}^{\alpha +1}+{\left(n-j\right)}^{\alpha +1},j=2,\dots n,\\ 1,j=n+1.\end{array}\right.$ (4.10)

• Explicit Direct Method using Quadratic Spline (EDMQS):

The Caputo fractional derivative can also be approximated by segmentary interpolation. Let [a, b] be an interval and $h=\frac{b-a}{n}$, we consider n + 1 nodes t _k = a + h.k, for k = 0, ..., n. Let $y:\left[a,b\right]\to {\mathrm{ℝ}}^m$ be a continuous function and y _k = y(t _k ), k = 0, ..., n. In Ferrari et al. ^9,11 the quadratic segmentary interpolator S(t) ∈ ℝ^m was determined such that it interpolates y(t) in t _k , k = 0, ..., n, S′(t) is continuous over the nodes t _k , k = 1, ..., n − 1 and also the quadratic deviation between the linear segmentary interpolation and the spline in the interval (t ₀ , t _n ) is minimum. S(t) is the function defined as: for k = 1, ..., n and t ∈ [t _k −1, t _k ]

$S\left(t\right)=\frac{t-t_{k-1}}{h}{y}_k-\frac{t-t_k}{h}{y}_{k-1}+a_k\left(t-t_{k-1}\right)\left(t-t_k\right)$ (4.11)

The coefficients a _k ∈ ℝ^m are written as $a_k=\sum _{j=0}^n{c}_{k,j}{y}_j$, where c _k,j is defined by ^9,11:

$c_{k,j}=\left\{\begin{array}{l}\frac{{\left(-1\right)}^k}{h^2}{\beta }_{1}ifj=0\\ \frac{{\left(-1\right)}^{k+1}}{h^2}\left(2{\beta }_1+{\beta }_2\right)ifj=1\\ \frac{{\left(-1\right)}^{k+j}}{h^2}\left({\beta }_{j-1}+2{\beta }_j+{\beta }_{j+1}\right)if1<j<n-1\\ \frac{{\left(-1\right)}^{k+n-1}}{h^2}\left(2{\beta }_{n-2}+{\beta }_{n-1}\right)ifj=n-1\\ \frac{{\left(-1\right)}^{k+n}}{h^2}{\beta }_{n-1}ifj=n\end{array}\right.$ (4.12)

with β _j = j/n if j ≤ k− 1 and β _j = j/n− 1 if j > k− 1. Here we set initial time a = 0 and m = 4.

Considering the definition of the Caputo derivative, the following approximation is obtained in Ferrari et al. ^8,10:

$D_C^{\alpha }y\left(t_k\right)\approx \frac{h^{1-\alpha }}{\Gamma \left(3-\alpha \right)}\sum _{j=1}^k{\gamma }_{k,j},$ (4.13)

being:

$$\begin{gathered} {\gamma }_{k,j}=2h\left[{\left(k-j\right)}^{1-\alpha }\left(j\alpha -j-k\right)+{\left(1-j+k\right)}^{1-\alpha }\left(-1+j+k+\alpha -\alpha j\right)\right]a_j \\ +\left(\alpha -2\right)\left[{\left(k-j\right)}^{1-\alpha }-{\left(1-j+k\right)}^{1-\alpha }\right]\left[h\left(1-2j\right)a_j+\frac{x_j-x_{j-1}}{h}\right]. \end{gathered}$$

From the above, given $\overline{Y}\left(t\right)$ an explicit and simple way to evaluate the Caputo derivative at the nodes t _k of the interval [0, t _final ] is obtained.

Taking into account that $\overline{Y}$ is approximated by S, (4.1) is evaluated in the nodes t _k , k = 1, ..., n which establishes a square algebraic system of order 4n in ${\left\{Y_{j,k}\right\}}_{j=1,k=1}^{4,n}$, with $Y_{j,k}={\overline{Y}}_j\left(t_k\right)$. Once the algebraic system is solved, the coefficients a _k of the spline S are determined by (4.12) and in this way the segmentary approximation S(t) is obtained.

4.1 Convergence results

In Ferrari et al. ¹⁴ the numerical methods GEEM (4.2), GIEM (4.3) and GTM (4.4) were implemented for different values of ∆t, considering fractional derivative order α = 0.99, reverse transcriptase inhibitor efficacy η = 0.6 and initial condition ${\overline{Y}}_j\left(0\right)=\left(300,10,10,10\right)$ for a total time of 120 days. The values and units of measurement of the parameters corresponding to the biology of the problem proposed in (3.1) are the usual ones in the biological literature and can be seen in Table 1.

Table 1:
Values of the parameters.

Parameter	Value	Unit
s	10	mm −3 day −1
k	0.000024	mm 3 day −1
µ	0.01	day −1
ε	0.4	day −1
b	0.05	day −1
µ 1	0.015	day −1
δ	0.26	day −1
N	1000	dimensionless
c	2.4	day
T (0)	300	mm −3
I(0)	10	mm −3
V (0)	10	mm −3
L(0)	10	mm −3
α	0.99	dimensionless

Figure 1 shows the graphs obtained by using these numerical methods corresponding to the ${\overline{Y}}_4$ component of the approximate solution, that is, to the variable L (virus density) as a function of time for each of the described schemes.

Figure 1:
Convergence study for GEEM, GTM and GIEM.

In Table 2, where the approximate solution obtained by each of the mentioned methods is indicated with ${\overline{Y}}_{\Delta t}$ taking ∆t as time step, the relative norms corresponding to each of the methods are detailed, calculated from the graphs just shown.

Table 2:
Relative norms for GEEM, GTM and GIEM.

	$a=\frac{||{\overline{Y}}_{{10}^{-1}}-{\overline{Y}}_{{10}^{-2}}||2}{||{\overline{Y}}_{{10}^{-1}}||2}$	$b=\frac{||{\overline{Y}}_{{10}^{-2}}-{\overline{Y}}_{{10}^{-3}}||2}{||{\overline{Y}}_{{10}^{-2}}||2}$	$\frac{a}{b}$
GEEM	0.017947345	0.017687582	0.9855264
GTM	0.000073027	0.000000621	0.0085037
GIEM	0.001921246	0.000195635	0.1018271

It is easy to see that the GTM and GIEM methods converge much faster than the GEEM method. The graphs obtained by GTM and GIEM are very similar and both schemes require practically the same implementation time. However, the convergence of GTM is better than that of GIEM and also the latter requires α ∈ (0, 1) while the former additionally allows α = 1.

4.2 Sensitivity results

4.2.1 Sensitivity results with respect to derivative order

In Ferrari et al. ¹⁴, we studied the sensitivity of the viral load L with respect to the fractional derivative order (Figure 2) using the GTM and η = 0. 6. Values of the fractional derivative order between α = 0. 5 and α = 1 were considered for a period of one year. It is clearly seen that L is not monotonic with respect to α. However, with the chosen values, after approximately 8 months no difference is noticed among the different outputs for α ≥ 0. 7, since L stabilizes at approximately 2700 units per mm ³. For times less than 8 months it is observed that for α chosen between 0. 96 and 1, the values of L are quite similar. This is profitable in the sense that although it is not clear which is the fractional derivative order that best fits the model to the reality, whichever of them will produce good approximations. While the classical model was good, the consideration of a fractional model produces corrections of the idealizations made when formulating the model.

Figure 2:
Sensitivity of L with respect to α using GTM (η = 0.6).

In Ferrari et al. ¹⁵, the fractional order numerical schemes EFEM, IFEM, FWDM and AM are considered for given values of α in the interval (0, 1). These numerical methods are known to be stable and have a linear order of convergence ⁵. Then, considering some values of α (the fractional derivative order), the component L of the approximate solution of the problem posed in (3.1) is studied again. To do this, ∆t = 0. 3 is used for the time step, η = 0. 6 and n = 2000 is the number of temporary nodes, that is, the evolution of the virus density over 600 days of drug application. This allows us to have a better estimate of the action of a drug on the patient in the long term, as will be seen later. In Figures 3, 4, 5 and 6, it can be seen that all the mentioned schemes behave in the same way in the first stages and they coincide quickly. In consequence any of them can be chosen to carry out the planned studies.

Figure 3:
L with α = 0. 6 (AM).

Figure 4:
L with α = 0. 7 (AM).

Figure 5:
L with α = 0. 8 (AM).

Figure 6:
L with α = 0. 9 (AM).

Since the AM is a scheme that improves the EFEM and avoids solving a system of equations at each instant of time like IFEM and FWDM, this is the numerical scheme that will be used to establish how the fractional derivative order α affects the model set in (3.1). This can be seen in Figures 3, 4, 5, 6 and 7.

Figure 7:
L with α between 0. 9 and 1 (AM).

Finally, in Ferrari et al. ¹², the methods AM and EDMQS applied to (3.1) are considered. EDMQS is the last numerical scheme with which the behavior of the numerical solution is investigated with respect to the fractional derivative order α varying between 0. 6 and 0. 9 for a total time of 600 days. The parameters mentioned in (3.1) and the initial conditions take the values of the Table 1.

The Figures 8, 9, 10 and 11 show the graphs of the variable of greatest interest, the viral load L, obtained as a result of the implementation of the mentioned numerical methods (LA is the one obtained with AM and L _Spline is the one obtained with EDMQS). The number of temporary nodes used in each of the schemes is not the same. The AM approximates the fractional derivative using a trapezoidal formula in conjunction with an implicit finite difference scheme (in this case n = 2000 nodes were used). On the other hand, the EDMQS needs the resolution of an algebraic system of order 4n that requires a great computational cost, so few temporary nodes were used for it (n = 100 nodes).

Figure 8:
L with α = 0. 6 (EDMQS).

Figure 9:
L with α = 0. 7 (EDMQS).

Figure 10:
L with α = 0. 8 (EDMQS).

Figure 11:
L with α = 0. 9 (EDMQS).

In all cases, the similarity of the behavior of the approximate solutions can be noted, except for a small difference in the first days of treatment. The advantage of AM over EDMQS is that it requires a lower computational cost and therefore can be implemented with a larger number of temporary nodes. However, we observe that EDMQS, despite having a higher computational cost, with less temporary nodes achieves very good precision, which is not the case with AM.

4.2.2 Sensitivity results with respect to drug efficacy

In Ferrari et al. ¹⁴, considering the critical value of the coefficient that represents the reverse transcriptase inhibitor efficacy η _crit = 0.88375, the behavior of the viral load L was studied for drug efficacy values η close to the mentioned η _crit (Figure 12). In fact, it is observed that for values of η < η _crit , even though the virus density decreases to almost 0, it increases again in longer instants of time as η increases. On the other hand, for η = 0.9, L tends to 0. Furthermore, it is observed that the maximum viral load decreases as the efficacy of the reverse transcriptase inhibitor increases, which was to be expected according to the biological interpretation of the problem. However, L is not monotonous with respect to η.

Figure 12:
Sensitivity of L with respect to η using GTM (α = 0.99).

In Ferrari et al. ¹⁵, we studied how the approximate solution of the initial value problem modeled in (3.1) behaves when varying the drug efficacy value η, given a fractional derivative order α = 0. 7. In Figure 13 the graph of L for different values of η is shown. We observe that for η = 0. 6 and η = 0. 65 the viral load increases up to stabilize in the infectious status. For the values η = 0. 7, η = 0. 75, η = 0. 8 and η = 0. 85, although in the first days the viral load decreases, then increases again, this growth occurring in longer moments of time as η grows. In addition, as is biologically expected, it is again observed that the maximum viral load decreases as the efficacy of the drug increases. On the other hand, from ηcrit = 0. 88375 the virus density gets closer and closer to 0, which is to be expected since in these cases the non-infectious status is asymptotically stable.

Figure 13:
L with α = 0. 7, η between 0. 6 and 0. 9 using AM.

5 COMPARISON WITH REAL DATA

In Ferrari et al. ¹⁵, we considered medical records of 10 patients selected from a total of 144. Patients were selected taking into account that a large amount of data were available over time from the start of treatment (at least 4 data) and that the person had good adherence to treatment. In that work, we found the derivative order α to use so that the numerical results obtained for the model described in (3.1) best fit those data. So, this determines whether a fractional derivative order or a classical derivative order is better.

For the classical model, i.e. the model (3.1) with α = 1, a numerical scheme described as follows was implemented: Considering the ordinary differential equation

$\overline{Y}\text{'}=\overline{f}\left(t,\overline{Y}\right),$ (5.1)

it is integrated into the time interval [t _n ,t _n+1 ]:

${\overline{\mathbf{Y}}}_{n+1}-{\overline{\mathbf{Y}}}_n={\int }_{t_n}^{t_{n+1}}\overline{f}\left(s,\overline{Y}\left(s\right)\right)ds.$ (5.2)

Next, by using a linear interpolation function for $\overline{f}$ between the coordinate points $\left(t_n,\overline{\mathbf{Y}}\left(t_n\right)\right)$ and $\left(t_{n+1},\overline{\mathbf{Y}}\left(t_{n+1}\right)\right)$, an implicit scheme is obtained:

${\overline{\mathbf{Y}}}_{n+1}-\frac{\Delta t}{2}\overline{f}\left(t_{n+1},{\overline{\mathbf{Y}}}_{n+1}\right)={\overline{\mathbf{Y}}}_n+\frac{\Delta t}{2}\overline{f}\left(t_n,{\overline{\mathbf{Y}}}_n\right).$ (5.3)

Now, we consider the Adams method (4.8) for values of the fractional derivative order α = 0. 5, 0. 6, 0. 7, 0. 8, 0. 9 and the scheme defined in (5.3) for the classical order, that is α = 1. In all cases, the first data from the patient in the treatment period was used as the initial condition $\left(T\left(0\right),I\left(0\right),V\left(0\right),L\left(0\right)\right)$ with a time step ∆t = 0. 25 (6 hs.). The efficacy considered was η = 0.936, since it corresponds to the efficacy of 3 drugs combined (as usual in most treatments), each with an efficacy of 0.6.

In the Figures 14, 15, 16 and 17 we show the graphs corresponding to the data of four of the patients considered, in which we compare the CD4⁺ (i.e. T +I+V ) of their medical records with the approximate solutions of (3.1) produced by the method for different values of α.

Figure 14:
Evolution of CD4⁺ for patient 1.

Figure 15:
Evolution of CD4⁺ for patient 3.

Figure 16:
Evolution of CD4⁺ for patient 6.

Figure 17:
Evolution of CD4⁺ for patient 9.

Different behaviors are observed according to the patient: each of them has a very particular evolution of their CD4⁺ with respect to the rest. In Figure 14 it is observed that during the first 800 days patient 1 evolved with a slope similar to the slope obtained for α = 0. 7. In Figure 15 it can be seen that patient 3 had an evolution that is maintained between the solutions obtained for α = 0. 7 and α = 0. 8. In Figure 16 it can be seen that α = 0. 6 is a good order of derivation to approximate the evolution of patient 6. Finally, in Figure 17 it is seen that in the first 250 days, patient 9 had an evolution that exceeds any of the approximate ones, but then the behavior is quite similar to that of α = 0. 9.

All this is reflected in Table 3 in which the averages of the relative errors $e_r=\frac{\left|CD{4}_{patient}-CD{4}_{approximated}\right|}{CD{4}_{patient}}$ are compared for each patient and each value of α. The smallest relative error in each row is highlighted in green.

Table 3:
Relative errors comparison.

Patient	α = 0. 5	α = 0. 6	α = 0. 7	α = 0. 8	α = 0. 9	α = 1
1	0. 3116	0. 1734	0. 1537	0. 2624	0. 4089	0. 4896
2	0. 136	0. 1051	0. 0674	0. 0542	0. 0629	0. 0828
3	0. 2503	0. 1712	0. 0742	0. 0652	0. 1782	0. 2623
4	0. 3626	0. 2925	0. 2022	0. 1098	0. 0535	0. 0515
5	0. 2611	0. 2307	0, 1908	0. 1503	0. 1261	0. 1327
6	0, 2091	0. 1858	0. 3875	0. 691	0. 9741	1. 1737
7	0. 1482	0. 1163	0. 0771	0. 0967	0. 1159	0. 1343
8	0. 0908	0. 0899	0. 141	0. 2346	0. 3101	0. 3591
9	0. 2283	0. 203	0. 1681	0. 1288	0. 0966	0. 1061
10	0. 1633	0. 1659	0. 1842	0. 2096	0. 2379	0. 2652
Average	0. 2161	0. 1734	0. 1646	0. 2003	0. 2564	0. 3057

In the last row of Table 3 the relative errors were averaged for each value of α. We see that among the derivative orders analyzed, the one that best fits the solution obtained by the numerical methods to the patient data in average behavior is α = 0. 7.

6 CONCLUSIONS

In Arafa et al. ³ a fractional order mathematical model was presented to study the effects of medicinal treatment in people infected with the human immunodeficiency virus (HIV). This model evaluates, among other magnitudes, the density of healthy and HIV-infected CD4⁺ T cells.

These data are very necessary for the subject infected by the virus given the effects that an antiretroviral treatment causes in it.

The objective of this work was to consider several numerical schemes that involved fractional derivative in order to compare their behaviors and obtain a good approximation of the mentioned model solution.

Comparing the numerical solutions obtained with GIEM, GEEM and GTM, for several time steps, we noted that GTM and GIEM converge much faster than GEEM. The graphs obtained by GTM and GIEM were very similar and both schemes required practically the same implementation time. However, the convergence of GTM was better than that of GIEM and also the latter required α ∈ (0, 1) while the former additionally allowed α = 1.

Also, the behavior of the approximate solutions was analyzed by considering EFEM, IEFM, FWDM, AM and EDMQS described in Section 4 for different fractional derivative orders. It was obtained that these numerical solutions tend to the approximate numerical solution of the classical problem (i.e., α = 1) when the fractional derivative order tends to 1. This was done in two stages. First AM was used, since it is a scheme that improves EFEM and avoids solving a system of equations at each instant of time like IEFM and FWDM. Then we also solved the problem using EDMQS and compared this numerical solution with AM. We concluded that AM can be implemented with a large number of temporary nodes and low computational cost but EDMQS with few temporary nodes already achieves a very good accuracy.

Following, the sensitivity of the L component (virus density) with respect to the drug efficacy parameter η (drug efficacy) was studied for values of η close to the critical value η _crit for the case α = 0.99 with GTM and for the case α = 0.7 with AM. In both cases, as is biologically expected, it is observed that the maximum viral load decreases as the efficacy of the drug increases. In addition, for η > η _crit = 0.88375 the viral load tends to 0, which is to be expected since in these cases the non-infectious status is asymptotically stable.

Finally, we studied which derivative order provided a numerical solution to the model (3.1) that is closest to what is observable in some clinical histories collected from patients. It was concluded that the best performing numerical solution is obtained for α = 0. 7. This confirms that a fractional derivative order is better than a classical order to study this type of biological problems, where memory processes are involved.

Acknowledgments

This work was partially supported by Universidad Nacional de Rosario through the project 80020180100064UR “Esquemas numéricos para la resolución de modelos de orden fraccionario relativos al tratamiento de la infección por VIH”. The first author was also supported by CONICET through a PhD fellowship.

REFERENCES

• ¹
  T. Anastasio. The fractional-order dynamics of brainstem vestibulo-oculomotor neurons. Byol. Cybernet., 72(1) (1994), 69-79.
• ²
  A. Arafa, S. Rida & M. Khalil. The effect of anti-viral drug treatment of human immunodeficiency virus type 1 (HIV-1) described by a fractional order model. Appl. Math. Model., 37(4) (2013), 2189-2196.
• ³
  A. Arafa, S. Rida & M. Khalil. A fractional-order model of HIV infection with drug therapy effect. J. Egyptian Math. Soc., 22(3) (2014), 538-543.
• ⁴
  D. Baleanu, K. Diethelm, E. Scalas & J. Trujillo. “Fractional Calculus: Models and Numerical Methods”, volume 3. World Scientific, London (2012).
• ⁵
  L. Changpin & Z. Fanhai. “Numerical Methods for Fractional Calculus”. Taylor and Francis Group (2015).
• ⁶
  K. Diethelm. “The Analysis of Fractional Differential Equations. An Application-Oriented Exposition Using Differential Operators of Caputo Type”. Springer-Verlag (2004).
• ⁷
  V. Djordjević, J. Jarić, B. Fabry, J. Fredberg & D. Stamenović. Fractional derivatives embody essential features of cell rheological behavior. Ann. Biomed. Eng., 31(6) (2003), 692-699.
• ⁸
  A. Ferrari, M. Gadella, L. Lara & E.S. Marcus. Approximate solutions of one dimensional systems with fractional derivative. International Journal of Modern Physics C, 31(7) (2020).
• ⁹
  A. Ferrari, L. Lara & E.S. Marcus. Interpolación por splines cuadráticos: obtención de una fórmula explícita. Mecánica Computacional (MECOM 2018), (2018), 1099-1108.
• ¹⁰
  A. Ferrari, L. Lara & E.S. Marcus. Resolución de ecuaciones fraccionarias mediante spline cuadrático. Proceedings of the VII MACI 2019, (2019), 197-200.
• ¹¹
  A. Ferrari, L. Lara & E.S. Marcus. Convergence analysis and parity conservation of a new form of a quadratic explicit spline with applications to integral equations. Journal of the Egyptian Mathematical Society, 28(30) (2020).
• ¹²
  A. Ferrari, L. Lara, M. Olguin & E.S. Marcus. Aplicación de dos métodos numéricos a un modelo de orden fraccionario para el tratamiento del VIH. Proceedings of the VIII MACI 2021, (2021), 55-58.
• ¹³
  A. Ferrari & E.S. Marcus. Study of a fractional-order model for HIV infection of CD4+ T-cells with treatment. Journal of Fractional Calculus and Applications, 11(2) (2020), 12-22.
• ¹⁴
  A. Ferrari, M. Olguin & E.S. Marcus. Resultados Numéricos para un Modelo de Orden Fraccionario del Tratamiento de la Infección por VIH. Proceedings of the VI MACI 2017, (2017), 9-12.
• ¹⁵
  A. Ferrari, M. Olguin & E.S. Marcus. Estudio comparativo de algunos esquemas numéricos para un modelo de orden fraccionario del tratamiento de la infección por VIH. Mecánica Computacional (MECOM 2018), (2018), 1749-1758.
• ¹⁶
  T. Hartley, C. Lorenzo & H.K. Qammer. Chaos in fractional order Chua’s system. IEEE Trans. Circuit Syst. I, 42(8) (1995), 485-490.
• ¹⁷
  N. Heymans. Fractional Calculus Description of Non-Linear Viscoelastic Behaviour of Polymers. Nonlinear Dynam., 38(1) (2004), 221-231.
• ¹⁸
  C. Li & F. Zeng. “Numerical Methods for Fractional Calculus”. Taylor and Francis Group, London (2015).
• ¹⁹
  Z. Odibat & S. Momani. An algorithm for the numerical solution of differential equations of fractional order. Journal of Applied Mathematics and Informatics, 26 (2008), 15-27.
• ²⁰
  K. Oldham. A signal-independent electroanalytical method. Anal. Chem., 44(1) (1972), 196-198.
• ²¹
  T. Poinot & J. Trigeassou. Identification of Fractional Systems Using an Output-Error Technique. Nonlinear Dynam., 38 (2004), 133-154.
• ²²
  F. Rihan. Numerical Modeling of Fractional-Order Biological Systems. Abstract and Applied Analysis, 2013 (2013).
• ²³
  P. Srivastava, M. Banerjee & P. Chandra. Modeling the drug therapy for HIV infection. Journal of Biological Systems, 17(2) (2009), 213-223.

Publication Dates

• Publication in this collection
  14 Nov 2022
• Date of issue
  Oct-Dec 2022

History

• Received
  30 Sept 2021
• Accepted
  24 Mar 2022