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

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: (1.1) dTdt=s−μTT+rT(1−T+ITmax)−k1VT, \frac{dT}{dt}=s-\mu_TT+rT\left(1-\frac{T+I}{T_{max}}\right)-k_1VT, (1.2) dIdt=k1′VT−μII, \frac{dI}{dt}=k_1'VT-\mu_II, (1.3) dVdt=NμbI−k1VT−μVV, \frac{dV}{dt}=N\mu_bI-k_1VT-\mu_VV, subject to conditions (1.4) 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 T_max is their carry capacity. The parameter k₁ represents the rate of infection of T-cells with free virus, and is the source term for infected cells. k1′ k_1' is the rate at which infected cells become actively infected, μ₁ 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.

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

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,t²,...} be a complete set of base functions and define V, the solution space, by taking all possible linear combinations of the set S_∞ as (2.1) 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, (2.2) 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} (2.3) 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 S_R = {ψ₀(t), ψ₁(t),...} be an infinite set of base functions that are linearly independent, and define (2.4) 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 S_R.

Now we define three nonlinear operators 𝒩₁,𝒩₂,𝒩₃ : V → U as (2.5) N1[T(t),I(t),V(t)]=dTdt−s+μTT−rT(1−T+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, (2.6) N2[T(t),I(t),V(t)]=dIdt−k1′VT+μII, \mathcal{N}_2[T(t), I(t), V(t)]= \frac{dI}{dt}-k_1'VT+\mu_II, (2.7) N3[T(t),I(t),V(t)]=dVdt−Nμ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 Φ₁(t;q) from T₀(t) to T (t), Φ₂(t;q) from I₀(t) to I(t), and Φ₃(t;q) from V₀(t) to V (t) we construct so-called zeroth-order deformation equations (2.8) 0=(1−q)ℒ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)], (2.9) 0=(1−q)ℒ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)], (2.10) 0=(1−q)ℒ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 ℒ₁,ℒ₂,ℒ₃ are linear operators, T₀(t),I₀(t),V₀(t) ∈ V^* are initial guesses, c₀ 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 Φ₁(t;q),Φ₂(t;q) and Φ₃(t;q) with respect to q as (2.11) Φ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.12) Φ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, (2.13) Φ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, (2.14) 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), (2.15) 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 k^th-order homotopy-derivative operator (see [9] for details).

Now applying k^th-order homotopy-derivative to (2.8)–(2.9), we obtain k^th-order deformation equations of OHAM (2.16) ℒ1[Tk(t)−χkTk−1(t)]=c0δk−1ξ(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.17) ℒ2[Ik(t)−χkIk−1(t)]=c0δk−1ξ(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, (2.18) ℒ3[Vk(t)−χkVk−1(t)]=c0δk−1ξ(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 (2.19) χk={0, k≤1,1, k>1. \chi_k = \begin{cases}0, \quad k \leq 1, \\1, \quad k > 1.\end{cases} and (2.20) δ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 ℒ₁, ℒ₂ and ℒ₃ so that the corresponding form of equations (2.16)–(2.18) can be determined. Then the number of terms desired in T_k, I_k and V_k 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 (2.21) Tk(t)=χkTk−1+c0J[δk−11(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, (2.22) Ik(t)=χkIk−1+c0J[δk−12(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, (2.23) Vk(t)=χkVk−1+c0J[δk−13(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 (2.24) 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 𝒥 : U → V by (2.25) 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.

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 L²-norm of residual error functions we define square residual error functions as (3.1) Eξ[h,A]=∫01(Nξ[T^(t),I^(t),V^(t)])2dt for ξ=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 (3.2) Eξ[h,A]=1M+1∑j=0M(Nξ[T^(jM+1),I^(jM+1),V^(jM+1)])2 for ξ=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 (3.3) 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 T₀ = 1000, I₀ = 0, V₀ = 0.001, r = 0.03, μ_T = 0.02, μ_b = 0.24, μ_V = 2.4, k₁ = 2.4 × 10⁻⁵, k1'=2×10−5 k_1^{'} = 2 \times 10^{-5} , s = 10, T_max = 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

Fig. 2

Fig. 3

Fig. 4

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 T₀ = 1000, I₀ = 0, V₀ = 0.001, r = 0.03, μ_T = 0.02, μ_I = 0.26, μ_b = 0.24, μ_V = 2.4, k₁ = 2.4 × 10⁻⁵, k1'=2×10−5 k_1^{'} = 2 \times 10^{-5} , N = 500, s = 10, and T_max = 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

Fig. 6

Fig. 7

Fig. 8

Fig. 9

Fig. 10

Fig. 11

Fig. 12

Fig. 13

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⁻¹⁰ (see the Table 2 and Figure 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.