1. bookAHEAD OF PRINT
Journal Details
License
Format
Journal
First Published
01 Jan 2016
Publication timeframe
2 times per year
Languages
English
access type Open Access

A Method of Directly Defining the inverse Mapping for a HIV infection of CD4+ T-cells model

Published Online: 15 Oct 2020
Page range: -
Received: 07 Aug 2019
Accepted: 06 Oct 2019
Journal Details
License
Format
Journal
First Published
01 Jan 2016
Publication timeframe
2 times per year
Languages
English
Abstract

In 2015, Shijun Liao introduced a new method of directly defining the inverse mapping (MDDiM) to approximate analytically a nonlinear differential equation. This method, based on the Homotopy Analysis Method (HAM) was proposed to reduce the time it takes in solving a nonlinear equation. Very recently, Dewasurendra, Baxter and Vajravelu (Applied Mathematics and Computation 339 (2018) 758–767) extended the method to a system of two nonlinear differential equations. In this paper, we extend it further to obtain the solution to a system of three nonlinear differential equations describing the HIV infection of CD4+ T-cells. In addition, we analyzed the advantages of MDDiM over HAM, in obtaining the numerical results. From these results, we noticed that the infected CD4+ T-cell density increases with the number of virions N; but decreases with the blanket death rate μI.

Keywords

MSC 2010

Introduction

Mathematical models become an important tool in analyzing the dynamics of human immunodeficiency virus (HIV), hepatitis B virus (HBV), and hepatitis C virus (HCV) infection [1]–[2]. The HIV mainly targets a host's CD4+ T-cells, which are the most abundant white blood cells of the immune system. Further, HIV wreaks the most havoc on the CD4+ T-cells by causing their destruction and decline, and decreasing the body ability to fight infection.

Chronic HIV infection causes gradual depletion of the CD4+ T-cell pool. Hence, progressively compromises the host's immune response to opportunistic infections, leading to acquired immunodeficiency syndrome (AIDS). For this reason, the count of CD4+ T-cells is a primary indicator used to measure progression of HIV infection.

In 1989, Perelson et al. [3] developed a simple model for the infection between the human immune system and HIV. Perelson et al. [4] then extended further the model, and observed that the model exhibits many of the symptoms of AIDS seen clinically: The long latency period, low levels of free virus in the body, and the depletion of CD4+ T-cells. They defined the model by considering four separate variables: Uninfected cells, latently infected cells, actively infected cells, and the concentration of free virus particles. Then, described the dynamics of these populations by a system of four coupled nonlinear differential equations. This model was then simplified by Culshaw and Ruan [5] by assuming that all of the infected cells are capable of producing the virus. In other words, they combined the latently infected cells and the actively infected ones.

This simplification resulted in a new system of three differential equations, which we will approximate using the expanded form of MDDiM. The simplified model is as follows: dTdt=sμTT+rT(1T+ITmax)k1VT,\frac{dT}{dt}=s-\mu_TT+rT\left(1-\frac{T+I}{T_{max}}\right)-k_1VT,dIdt=k1VTμII,\frac{dI}{dt}=k_1'VT-\mu_II,dVdt=NμbIk1VTμVV,\frac{dV}{dt}=N\mu_bI-k_1VT-\mu_VV, subject to conditions T(0)=T0,I(0)=I0,V(0)=V0.T(0)=T_0,\quad I(0)=I_0,\quad V(0)=V_0.

In this model, the variables T(t), I(t), and V(t) are the uninfected CD4+ T-cells, the infected CD4+ T-cells, and the concentration free HIV particles at time t respectively. Here, s is the source of the CD4+ T-cells from precursors, μT is the natural death rate of CD4+ T-cells, r is their growth rate (r > μT in general), and Tmax is their carry capacity. The parameter k1 represents the rate of infection of T-cells with free virus, and is the source term for infected cells. k1k_1' is the rate at which infected cells become actively infected, μ1 is a blanket death term for infected cells, μb is the lytic death rate for infected cells, μV is the loss rate of virus, N is the number of viral particles are released by each lysing cell. Table 1 summarizes parameters and variables.

Variables and parameters for viral spread

Parameters and VariablesMeaning
DependentVariables
TUninfected CD+ T-cell population size
IInfected CD+ T-cell density
VInitial density of HIV RNA
Parameters and Constants
μTNatural death rate of CD+ T-cell
μIBlanket death rate of infected CD+ T-cel
μbLytic death rate for infected cells
μvDeath rate of free virus
k1Rate CD+ T-cel become infected with virus
k1k_1'Rate infected cells becomes active
rGrowth rate of CD+ T-cell population
NNumber of virions produced by infected CD+ T-cel
TmaxMaximal population level of CD+ T-cel
sSource term for uninfected CD+ T-cel
Derived quantities
T0CD+ T-cel population for HIV-negative persons

This particular model and its variations have been solved multiple times by different methods including the HAM by M. Goreishi et al. [6], Homotopy Perturbation Method [7] and the Laplace Adomian decomposition method (LADM) by Ongun [8]. In the present paper we solve the model by MDDiM. This is the first time we used MDDiM to solve a system of three nonlinear coupled ordinary differential equations arise in the field of mathematical biology, using a single inverse linear map to solve the three deformation equations (2.21)–(2.23).

HAM and MDDiM approach

In this section, we first discuss the space that solution and base functions come from, and then will introduce the Method of Directly Defining inverse Mapping Method (MDDiM). MDDiM is an extension to the Optimal Homotopy Analysis Method (OHAM). In OHAM we solve an infinite number of linear ordinary differential equations to obtain series solution but in MDDiM with the help of directly defined inverse map 𝒥 we solve systems of linear equations.

Let S = {1,t,t2,...} be a complete set of base functions and define V, the solution space, by taking all possible linear combinations of the set S as V={k=0+aktk,ak}.V=\{ \sum_{k=0}^{+\infty}a_kt^k,\quad a_k\in \mathcal{R}\}. Next we define a space for the initial guesses by taking all possible linear combinations of the set S* = {1,t}. That is, V*={a0+a1t|ao,a1}.V^*=\{a_0+a_1 t\left | a_o, a_1\in \mathcal{R}\right. \}. Here we use the first two base functions of the set S to define S* because each equation in our coupled nonlinear system (1.1) has two initial conditions. Next, we define V^\hat{V}V^={k=2+aktk,ak}\hat{V}=\{\sum_{k=2}^{+\infty} a_kt^k, \quad a_k\in \mathcal{R} \} such that, V=V*V^V=V^* \cup \hat{V} .

Similarly, let SR = {ψ0(t), ψ1(t),...} be an infinite set of base functions that are linearly independent, and define U={m=1+bmψm(t)|bm}U=\{\sum_{m=1}^{+\infty}b_m\psi_m(t)\left| b_m\in\mathcal{R}\right.\} by taking all possible linear combinations of SR.

Now we define three nonlinear operators 𝒩1,𝒩2,𝒩3 : VU as N1[T(t),I(t),V(t)]=dTdts+μTTrT(1T+ITmax)+k1VT,\mathcal{N}_1[T(t), I(t), V(t)]=\frac{dT}{dt}-s+\mu_TT-rT\left(1-\frac{T+I}{T_{max}}\right)+k_1VT,N2[T(t),I(t),V(t)]=dIdtk1VT+μII,\mathcal{N}_2[T(t), I(t), V(t)]= \frac{dI}{dt}-k_1'VT+\mu_II,N3[T(t),I(t),V(t)]=dVdtNμb+k1VT+μVV.\mathcal{N}_3[T(t), I(t), V(t)]= \frac{dV}{dt}-N\mu_b+k_1VT+\mu_VV. Next, to build continuous variations Φ1(t;q) from T0(t) to T (t), Φ2(t;q) from I0(t) to I(t), and Φ3(t;q) from V0(t) to V (t) we construct so-called zeroth-order deformation equations 0=(1q)1[Φ1(t;q)T0(t)]c0qN1[Φ1(t;q),Φ2(t;q),Φ3(t;q)],0=(1-q)\mathcal{L}_1[\Phi_1(t;q)-T_0(t)]-c_0q\mathcal{N}_1[\Phi_1(t;q),\Phi_2(t;q),\Phi_3(t;q)],0=(1q)2[Φ2(t;q)I0(t)]c0qN2[Φ1(t;q),Φ2(t;q),Φ3(t;q)],0=(1-q)\mathcal{L}_2[\Phi_2(t;q)-I_0(t)]-c_0q\mathcal{N}_2[\Phi_1(t;q),\Phi_2(t;q),\Phi_3(t;q)],0=(1q)3[Φ3(t;q)V0(t)]c0qN3[Φ1(t;q),Φ2(t;q),Φ3(t;q)]].0=(1-q)\mathcal{L}_3[\Phi_3(t;q)-V_0(t)]-c_0q\mathcal{N}_3[\Phi_1(t;q),\Phi_2(t;q),\Phi_3(t;q)]]. Here 1,2,3 are linear operators, T0(t),I0(t),V0(t) ∈ V* are initial guesses, c0 is the convergence control parameter and q ∈ [0,1] is the homotopy parameter. The convergence control parameter will be used to optimize the function approximation in the next section.

Assuming that the solutions to the zeroth-order deformation equations are analytic at q = 0, we write Maclaurin series for Φ1(t;q),Φ2(t;q) and Φ3(t;q) with respect to q as Φ1(t;q)=T0(t)+k=1k=+Tk(t)qk,\Phi_1(t;q)=T_0(t)+\sum_{k=1}^{k=+\infty}T_k(t)q^k,Φ2(t;q)=I0(t)+k=1k=+Ik(t)qk,\Phi_2(t;q)=I_0(t)+\sum_{k=1}^{k=+\infty}I_k(t)q^k,Φ3(t;q)=V0(t)+k=1k=+Vk(t)qk,\Phi_3(t;q)=V_0(t)+\sum_{k=1}^{k=+\infty}V_k(t)q^k, where, Tk(t)=1k!kΦ1(t;q)qk|q=0=DkΦ1(q;t),Ik(t)=1k!kΦ2(t;q)qk|q=0=DkΦ2(q;t),T_k(t)=\left.\frac{1}{k!}\frac{\partial^k\Phi_1(t;q)}{\partial q^k}\right |_{q=0}=\mathscr{D}_k\Phi_1(q;t),\quad I_k(t)=\left.\frac{1}{k!}\frac{\partial^k\Phi_2(t;q)}{\partial q^k}\right |_{q=0}=\mathscr{D}_k\Phi_2(q;t),Vk(t)=1k!kΦ3(t;q)qk|q=0=DkΦ3(q;t).V_k(t)=\left.\frac{1}{k!}\frac{\partial^k\Phi_3(t;q)}{\partial q^k}\right |_{q=0}=\mathscr{D}_k\Phi_3(q;t). Here 𝒟k is called kth-order homotopy-derivative operator (see [9] for details).

Now applying kth-order homotopy-derivative to (2.8)–(2.9), we obtain kth-order deformation equations of OHAM 1[Tk(t)χkTk1(t)]=c0δk1ξ(t),Tk(0)=0,Tk(0)=0,\mathcal{L}_1[T_k(t)-\chi_k T_{k-1}(t)]=c_0\delta_{k-1}^{\xi}(t), \quad T_k(0)=0,\quad T_k'(0)=0, 2[Ik(t)χkIk1(t)]=c0δk1ξ(t),Ik(0)=0,Ik(0)=0,\mathcal{L}_2 [I_k(t)-\chi_k I_{k-1}(t)]=c_0\delta_{k-1}^{\xi}(t), \quad I_k(0)=0,\quad I_k'(0)=0,3[Vk(t)χkVk1(t)]=c0δk1ξ(t),Vk(0)=0,Vk(0)=0,\mathcal{L}_3[V_k(t)-\chi_k V_{k-1}(t)]=c_0\delta_{k-1}^{\xi}(t), \quad V_k(0)=0,\quad V_k'(0)=0, where χk={0,k1,1,k>1.\chi_k = \begin{cases}0, \quad k \leq 1, \\1, \quad k > 1.\end{cases} and δkξ(t)=Dk[Nξ[Φ1(t;q),Φ2(t;q),Φ3(t,q)]]\delta_{k}^{\xi}(t)= \mathscr{D}_k[\mathcal{N}_{\xi}[\Phi_1(t;q),\Phi_2(t;q),\Phi_3(t,q)]] for ξ = 1, 2, 3.

The benefit of OHAM is that we have a considerable freedom to choose linear operators 1, 2 and 3 so that the corresponding form of equations (2.16)–(2.18) can be determined. Then the number of terms desired in Tk, Ik and Vk of the series solutions can be determine iteratively. OHAM has been successfully used in science and engineering applications (see [10][11][12][13][14]), and the choice of axiliary linear operators and convergent control parameters were studied in [11]. The only drawback of OHAM is it takes lot of CPU time to calculate inverse linear operator. To overcome this obstacle, Liao [15] and Dewasurendra et al. [16] introduced the Method of Directly Defining inverse Mapping.

Now applying the inverse linear map to OHAM deformation equations (2.16)–(2.18). We obtain the following deformation equations in the frame of MDDiM Tk(t)=χkTk1+c0J[δk11(t)]+ak,0+ak,1t,T_k(t)=\chi_kT_{k-1}+c_0\mathscr{J}[\delta_{k-1}^{1}(t)]+a_{k,0}+a_{k,1}t,Ik(t)=χkIk1+c0J[δk12(t)]+bk,0+bk,1t,I_k(t)=\chi_kI_{k-1}+c_0\mathscr{J}[\delta_{k-1}^{2}(t)]+b_{k,0}+b_{k,1}t,Vk(t)=χkVk1+c0J[δk13(t)]+ck,0+ck,1t,V_k(t)=\chi_kV_{k-1}+c_0\mathscr{J}[\delta_{k-1}^{3}(t)]+c_{k,0}+c_{k,1}t, subject to the initial conditions Tk(0)=Ik(0)=Vk(0)=0,Tk(0)=Ik(0)=Vk(0)=0.T_k(0)=I_k(0)=V_k(0)=0, \quad T'_k(0)=I_k'(0)=V_k'(0)=0.

Here we use only one inverse linear operator to all three deformation equations which leads to less complicated solutions. However, a different inverse linear map could be used to obtain series solutions in different structures. Here we define the inverse linear map 𝒥 : UV by J[tk]=tk+1Ak+1,\mathscr{J}[t^k]=\frac{t^{k+1}}{Ak+1}, where A is a parameter which will be used to optimize the square residual error functions.

Results and Error Analysis

The approximate series solution for the coupled nonlinear system (1.1)–(1.3) with boundary conditions (1.4) are obtained using MDDiM. Further, error analysis is carried out to get a general idea about how accurate the approximate solutions are.

First, define three term approximations for T^\widehat{T} , I^\widehat{I} and V^\widehat{V} which are sum of the first four terms to the MDDiM deformation equations. Then, define residual error functions N1[T^(t),I^(t),V^(t)]\mathcal{N}_1[\widehat{T}(t),\widehat{I}(t),\widehat{V}(t)] , N2[T^(t),I^(t),V^(t)]\mathcal{N}_2[\widehat{T}(t),\widehat{I}(t),\widehat{V}(t)] and N3[T^(t),I^(t),V^(t)]\mathcal{N}_3[\widehat{T}(t),\widehat{I}(t),\widehat{V}(t)] in order to obtain square residual error functions. Now, taking the square of the L2-norm of residual error functions we define square residual error functions as Eξ[h,A]=01(Nξ[T^(t),I^(t),V^(t)])2dtforξ=1,2,3.E_{\xi}[h, A]=\int_0^{1}\left(\mathcal{N}_{\xi}\left[\widehat{T}(t),\widehat{I}(t),\widehat{V}(t)\right]\right)^2dt \quad\text{for}\quad \xi=1,2,3. However, in practice the evaluation of Eξ [h,A,β] is much time consuming so instead of exact residual error we use average residual error defined as Eξ[h,A]=1M+1j=0M(Nξ[T^(jM+1),I^(jM+1),V^(jM+1)])2forξ=1,2,3,E_{\xi}[h,A]=\frac{1}{M+1}\sum_{j=0}^M\left(N_{\xi}\left[\widehat{T}\left(\frac{j}{M+1}\right),\widehat{I}\left(\frac{j}{M+1}\right),\widehat{V}\left(\frac{j}{M+1}\right)\right]\right)^2\quad \text{for}\quad \xi=1,2,3, and define a total error function by taking affine combination of three square residual functions E[h,A]=ξ=13Eξ[h,A].E[h,A]=\sum_{\xi=1}^3 E_{\xi}[h,A].

Next, we optimize the total error function with respect to h and A and obtain optimal values for h and A. Substituting optimal values of h and A we obtain three term approximation solution to the nonlinear coupled system (1.1)–(1.3) which satisfy the conditions (1.4).

Taking three sets of parametric values for μI and N while keeping T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , s = 10, Tmax = 1500 fixed, we obtained the minimum value of the total error function E[h,A] along with suitable values for h and A, and presented in Table 1. The plots of the total error functions E[h,A] for three sets of representative parametric values are plotted in Figures 1, 2 and 3 for t ∈ [1,499]. Further, residual error function E[h,A] verses number of terms of the series approximation is presented in Figure 4 for two sets of parameters to guarantee the convergence of the series solution. From these curves, it is clear that the solutions obtained by MDDiM gives an analytical solution with high accuracy only with few iterations.

Fig. 1

Plot of E(h,A) for T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500.

Fig. 2

Plot of E(h,A) for T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 600, s = 10, Tmax = 1500.

Fig. 3

Plot of E(h,A) for T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500.

Fig. 4

Residual Error verses Terms of approximation for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where curve 1 has N = 500, and Curve 2 has N = 600.

Solution curves for uninfected CD+ T-cell population size T, Infected CD+ T-cell density I and Initial density of HIV RNA (V) are presented in Figures 5–13 along with numerical results we obtained by Runge-Kutta method. Further, 4-term MDDiM solutions, 6-term OHAM solutions proposed by Ghoreishi et al., and numerical results obtained by Runge-Kutta method for the parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, and Tmax = 1500 are presented in table 3 - table 4. From these tables, it is clear that the MDDiM gives an analytical solution with high accuracy only with few iterations when compare to OHAM solution.

Fig. 5

Comparison of the MDDiM solution of T (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 6

Comparison of the MDDiM solution of I(t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 7

Comparison of MDDiM solution of V (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 8

Comparison of MDDiM solution of T (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 600, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 9

Comparison of MDDiM solution of I(t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 600, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 10

Comparison of MDDiM solution of V (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 600, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 11

Comparison of MDDiM solution of T (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 12

Comparison of MDDiM solution of I(t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 13

Comparison of MDDiM solution of V (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Conclusions

In the present study, MDDiM has been developed and used to solve the model of HIV infection of CD4+ T-cells. Approximate analytical solutions for the uninfected CD4+T-cell population size, infected CD4+ T-cell density and initial density of HIV RNA were found. Our analytical solutions are in good agreement with the results of Ghoreishi et al. [4], and with the numerical results, we obtained by Runge-Kutta method (see Figures 5–13). From these results, we noticed that the infected CD4+ T-cell density increases with the number of virions N; but decreases with the blanket death rate μI (see Figures 6, 9 and 12).

Since, inverse mapping is directly defined, approximate series solutions are obtained with less CPU time compare with OHAM solution. Also, it is investigated that selected inverse map leads converge series solutions with total error 10−10 (see the Table 2 and Figure 4).

Minimum of the squared residual error E(h,A) for three different sets of μI and N for fixed parametric values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×105k_1^{'} = 2 \times 10^{-5} , s = 10, Tmax = 1500.

μINhAE[h,A]
0.26500−0.42950.48694.068 × 10−10
0.26600−0.41870.47315.682 × 10−10
0.1500−0.43110.48524.235 × 10−10

Numerical comparison for I(t)

tMDDiMOHAMRunge-Kutta
0.0000
0.23.1568 × 10−61.7626 × 10−73.1318 × 10−6
0.45.1445 × 10−62.9451 × 10−75.1352 × 10−6
0.66.5569 × 10−63.7344 × 10−76.5894 × 10−6
0.87.7445 × 10−64.2602 × 10−77.8026 × 10−6
1.08.8153 × 10−64.6108 × 10−78.9418 × 10−6

Numerical comparison for V (t)

tMDDiMOHAMRunge-Kutta
0.00.010.010.01
0.26.4618 × 10−46.1744 × 10−45.5042 × 10−4
0.44.7805 × 10−43.8431 × 10−44.8189 × 10−4
0.64.0834 × 10−42.4258 × 10−44.0993 × 10−4
0.83.8591 × 10−41.5659 × 10−43.9044 × 10−4
1.03.9576 × 10−41.0406 × 10−44.0056 × 10−4

In this study, we used a single inverse linear operator for all three deformation equations (2.21)–(2.23), but it is open to use two, or three inverse maps instead. Further, convergence of the series solutions may depend on the choice of inverse linear map and also with the defined solution space. So, it is worth to investigate properties of inverse linear maps in the frame of MDDiM.

To the best of our knowledge, this is the first time someone has used this method to solve HIV infection model. This novel method is more general and can be used to analyze non-linear systems of differential equations arising in science and engineering problems.

Fig. 1

Plot of E(h,A) for T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 500, s = 10, Tmax = 1500.
Plot of E(h,A) for T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500.

Fig. 2

Plot of E(h,A) for T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 600, s = 10, Tmax = 1500.
Plot of E(h,A) for T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 600, s = 10, Tmax = 1500.

Fig. 3

Plot of E(h,A) for T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 500, s = 10, Tmax = 1500.
Plot of E(h,A) for T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500.

Fig. 4

Residual Error verses Terms of approximation for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 500, s = 10, Tmax = 1500, where curve 1 has N = 500, and Curve 2 has N = 600.
Residual Error verses Terms of approximation for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where curve 1 has N = 500, and Curve 2 has N = 600.

Fig. 5

Comparison of the MDDiM solution of T (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.
Comparison of the MDDiM solution of T (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 6

Comparison of the MDDiM solution of I(t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.
Comparison of the MDDiM solution of I(t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 7

Comparison of MDDiM solution of V (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.
Comparison of MDDiM solution of V (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 8

Comparison of MDDiM solution of T (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 600, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.
Comparison of MDDiM solution of T (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 600, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 9

Comparison of MDDiM solution of I(t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 600, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.
Comparison of MDDiM solution of I(t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 600, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 10

Comparison of MDDiM solution of V (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 600, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.
Comparison of MDDiM solution of V (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.26, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 600, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 11

Comparison of MDDiM solution of T (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.
Comparison of MDDiM solution of T (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 12

Comparison of MDDiM solution of I(t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.
Comparison of MDDiM solution of I(t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Fig. 13

Comparison of MDDiM solution of V (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, 
k1'=2×10−5k_1^{'} = 2 \times 10^{-5}
, N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.
Comparison of MDDiM solution of V (t) for parameter values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μI = 0.1, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, Tmax = 1500, where Curve 1 is Runge-Kutta results, Curve 2 is MDDiM results.

Minimum of the squared residual error E(h,A) for three different sets of μI and N for fixed parametric values T0 = 1000, I0 = 0, V0 = 0.001, r = 0.03, μT = 0.02, μb = 0.24, μV = 2.4, k1 = 2.4 × 10−5, k1'=2×10−5k.1^{'} = 2 \times 10^{-5} , s = 10, Tmax = 1500.

μINhAE[h,A]
0.26500−0.42950.48694.068 × 10−10
0.26600−0.41870.47315.682 × 10−10
0.1500−0.43110.48524.235 × 10−10

Numerical comparison for I(t)

tMDDiMOHAMRunge-Kutta
0.0000
0.23.1568 × 10−61.7626 × 10−73.1318 × 10−6
0.45.1445 × 10−62.9451 × 10−75.1352 × 10−6
0.66.5569 × 10−63.7344 × 10−76.5894 × 10−6
0.87.7445 × 10−64.2602 × 10−77.8026 × 10−6
1.08.8153 × 10−64.6108 × 10−78.9418 × 10−6

Numerical comparison for V (t)

tMDDiMOHAMRunge-Kutta
0.00.010.010.01
0.26.4618 × 10−46.1744 × 10−45.5042 × 10−4
0.44.7805 × 10−43.8431 × 10−44.8189 × 10−4
0.64.0834 × 10−42.4258 × 10−44.0993 × 10−4
0.83.8591 × 10−41.5659 × 10−43.9044 × 10−4
1.03.9576 × 10−41.0406 × 10−44.0056 × 10−4

Variables and parameters for viral spread

Parameters and VariablesMeaning
DependentVariables
TUninfected CD+ T-cell population size
IInfected CD+ T-cell density
VInitial density of HIV RNA
Parameters and Constants
μTNatural death rate of CD+ T-cell
μIBlanket death rate of infected CD+ T-cel
μbLytic death rate for infected cells
μvDeath rate of free virus
k1Rate CD+ T-cel become infected with virus
k1k_1'Rate infected cells becomes active
rGrowth rate of CD+ T-cell population
NNumber of virions produced by infected CD+ T-cel
TmaxMaximal population level of CD+ T-cel
sSource term for uninfected CD+ T-cel
Derived quantities
T0CD+ T-cel population for HIV-negative persons

X. Wang, X. Song, Global stability and periodic solution of a model for HIV infection of CD4+ T-cells, Applied Mathematics and Computation 189 (2007) 1331–1340.WangX.SongX.Global stability and periodic solution of a model for HIV infection of CD4+ T-cellsApplied Mathematics and Computation189200713311340Search in Google Scholar

R.V. Culshaw, S. Ruan, A delay-differential equation model of HIV infection of CD4+ T-cells, Mathematical Biosciences 165 (2000) 27–39.CulshawR.V.RuanS.A delay-differential equation model of HIV infection of CD4+ T-cellsMathematical Biosciences16520002739Search in Google Scholar

A.S. Perelson, Modelling the interaction of the immune system with HIV, in: C. Castillo-Chavez (Ed.), Mathematical and Statistical Approaches to AIDS Epidemiology, Springer, Berlin, 1989, p. 350.PerelsonA.S.Modelling the interaction of the immune system with HIVin:Castillo-ChavezC.(Ed.),Mathematical and Statistical Approaches to AIDS EpidemiologySpringerBerlin1989350Search in Google Scholar

A.S. Perelson, D.E. Kirschner, R. De Boer, Dynamics of HIV infection of CD4+ T-cells, Mathematical Biosciences 114 (1993) 81.PerelsonA.S.KirschnerD.E.De BoerR.Dynamics of HIV infection of CD4+ T-cellsMathematical Biosciences114199381Search in Google Scholar

R.V. Culshaw, S. Ruan, A delay-differential equation model of HIV infection of CD4+ T-cells, Mathematical Biosciences 165 (2000) 27–39.CulshawR.V.RuanS.A delay-differential equation model of HIV infection of CD4+ T-cellsMathematical Biosciences16520002739Search in Google Scholar

M. Ghoreishi, A.I.B.Md. Ismail, A.K. Alomari, Application of the homotopy analysis method for solving a model for HIV infection of CD4+ T-cells, Mathematical and Computer Modelling (2011) 3007–3015.GhoreishiM.IsmailA.I.B.Md.AlomariA.K.Application of the homotopy analysis method for solving a model for HIV infection of CD4+ T-cellsMathematical and Computer Modelling201130073015Search in Google Scholar

M. Merdan, Homotopy perturbation method for solving a model for HIV infection of CD4+ T-cells, İstanbul Ticaret Üniversitesi Fen Bilimleri Dergisi Yil: 6 Sayi: 12 Güz 2007/2 s. pp. 39–52.MerdanM.Homotopy perturbation method for solving a model for HIV infection of CD4+ T-cellsİstanbul Ticaret Üniversitesi Fen Bilimleri DergisiYil: 6 Sayi: 12 Güz200723952Search in Google Scholar

M.Y. Ongun, The Laplace Adomian decomposition method for solving a model for HIV infection of CD4+ T cells, Mathematical and Computer Modelling 53 (2011) 597–603.OngunM.Y.The Laplace Adomian decomposition method for solving a model for HIV infection of CD4+ T cellsMathematical and Computer Modelling532011597603Search in Google Scholar

S. Liao, Notes on the homotopy analysis method: Some definitions and theorems. Commun. Nonlinear Sci. Numer. Simul. 14(4), 983–997 (2009).LiaoS.Notes on the homotopy analysis method: Some definitions and theoremsCommun. Nonlinear Sci. Numer. Simul.1449839972009Search in Google Scholar

S. Liao, Beyond Perturbation: Introduction to the Homotopy Analysis Method, Chapman & Hall/CRC Press, Boca Raton, 2003.LiaoS.Beyond Perturbation: Introduction to the Homotopy Analysis MethodChapman & Hall/CRC PressBoca Raton2003Search in Google Scholar

M. Baxter, R.A. Van Gorder and K. Vajravelu, On the choice of auxiliary linear operator in the optimal homotopy analysis of the Cahn-Hilliard initial value problem, Numerical Algorithms 66 (2014) 269–298.BaxterM.Van GorderR.A.VajraveluK.On the choice of auxiliary linear operator in the optimal homotopy analysis of the Cahn-Hilliard initial value problemNumerical Algorithms662014269298Search in Google Scholar

Y. Tan and S. Abbasbandy, Homotopy analysis method for quadratic Riccati differential equation, Communications in Nonlinear Science and Numerical Simulation 13 (2008) 539–546.TanY.AbbasbandyS.Homotopy analysis method for quadratic Riccati differential equationCommunications in Nonlinear Science and Numerical Simulation132008539546Search in Google Scholar

S. Liao, An optimal homotopy-analysis approach for strongly nonlinear differential equations, Communications in Nonlinear Science and Numerical Simulation 15 (2010) 2315–2332.LiaoS.An optimal homotopy-analysis approach for strongly nonlinear differential equationsCommunications in Nonlinear Science and Numerical Simulation15201023152332Search in Google Scholar

M. Baxter, R.A. Van Gorder and K. Vajravelu, Optimal analytic method for the nonlinear Hasegawa-Mima equation, The European Physical Journal Plus 129 (2014) 98. https://doi.org/10.1140/epjp/i2014-14098-xBaxterM.Van GorderR.A.VajraveluK.Optimal analytic method for the nonlinear Hasegawa-Mima equationThe European Physical Journal Plus129201498https://doi.org/10.1140/epjp/i2014-14098-xSearch in Google Scholar

S. Liao, Y Zhao, On the method of directly defining inverse mapping for nonlinear differential equations, Numerical Algorithms 72(4) (2016) 989–1020.LiaoS.ZhaoYOn the method of directly defining inverse mapping for nonlinear differential equationsNumerical Algorithms72420169891020Search in Google Scholar

M. Dewasurendra, M. Baxter, K. Vajravelu, A method of directly defining the inverse mapping for solutions of nonlinear coupled systems arising in convection heat transfer in a second grade fluid, Applied Mathematics and Computation 339 (2018) 758–767.DewasurendraM.BaxterM.VajraveluK.A method of directly defining the inverse mapping for solutions of nonlinear coupled systems arising in convection heat transfer in a second grade fluidApplied Mathematics and Computation3392018758767Search in Google Scholar

Recommended articles from Trend MD

Plan your remote conference with Sciendo