Arshad et al. Advances in Difference Equations (2017) 2017:92 DOI 10.1186/s13662-017-1143-0

0 Advances in Difference Equations

a SpringerOpen Journal

RESEARCH

Open Access

Effects of HIV infection on CD4+ T-cell population based on a fractional-order model

Sadia Arshad1,6*©, Dumitru Baleanu2,4, Weiping Bu3 and Yifa Tang1,5

CrossMark

"Correspondence: sadia_735@yahoo.com 1LSEC, ICMSEC, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, 100190, China 6COMSATS Institute of Information Technology, Lahore, Pakistan Full list of author information is available at the end of the article

Abstract

In this paper, we study the HIV infection model based on fractional derivative with particular focus on the degree of T-cell depletion that can be caused by viral cytopathicity. The arbitrary order of the fractional derivatives gives an additional degree of freedom to fit more realistic levels of CD4+ cell depletion seen in many AIDS patients. We propose an implicit numerical scheme for the fractional-order HIV model using a finite difference approximation of the Caputo derivative. The fractional system has two equilibrium points, namely the uninfected equilibrium point and the infected equilibrium point. We investigate the stability of both equilibrium points. Further we examine the dynamical behavior of the system by finding a bifurcation point based on the viral death rate and the number of new virions produced by infected CD4+ T-cells to investigate the influence of the fractional derivative on the HIV dynamics. Finally numerical simulations are carried out to illustrate the analytical results.

Keywords: fractional derivative; HIV model; finite difference scheme; dynamical analysis

ft Spri

ringer

1 Introduction

According to WHO there were approximately 35 million people at the end of 2013 living with HIV with 2.1 million people becoming newly infected in 2013 globally with HIV. HIV belongs to the family of lentiviruses, which means being acting slowly. Lentiviruses cause diseases that progress over a long period disturbing the immune system in humans. HIV produces virus particles by converting viral RNA into DNA in the cell and then making many RNA copies. The transformation is completed with an enzyme named reverse transcriptase. The change from RNA to DNA and back to RNA is substantial and makes fighting against HIV challenging. There is a chance of the virus mutating and there being errors each time when it happens. These copies or virus particles destroy the cell after formation and infecting other cells. Although HIV attacks many different cells, it inflicts the most disorder on the CD4+ T-cells, the main target of CD4+ T-cells is to form the body's overall immune response to external infections. HIV provides the basis to grow acquired immune deficiency syndrome (AIDS). For a person infected from HIV it can take 10-15 years to develop AIDS. On the medical frontier there have been many advances, but still there is no effective cure or vaccine available for HIV. Since the early 1980s, many mathematical models have been developed to better understand the interaction of HIV and the human immune system for the purpose of testing treatment strate© The Author(s) 2017. This article is distributed under the terms of the Creative Commons Attribution 4.0 InternationalLicense (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use,distribution, and reproduction in any medium, provided you give appropriate credit to the originalauthor(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

gies [1-8]. Silva and Torres [9] proposed a population model for TB-HIV/AIDS coinfec-tion transmission dynamics, they applied optimal control theory to TB-HIV/AIDS model to study optimal strategies for the minimization of the number of individuals with TB and AIDS. Recently Rocha et al. [10] investigated an HIV optimal control problem with delays in both state and control variables, where the objective is to find the optimal treatment strategy that maximizes the number of CD4+ T-cells and CTL immune response cells, keeping the drug strength low. Luo et al. [11] studied bifurcations and stability of an HIV model that incorporates the immune responses. An HIV model including latent infection and antiretroviral therapy is examined by Wang et al. [12]. They obtained the global asymptotic stability of the uninfected equilibrium by constructing a Lyapunov function. We refer the reader to the excellent review paper on mathematical modeling of HIV on different phenomena of [13]. Quantitative analysis of HIV-1 replication in vivo has made significant contributions to understanding of AIDS pathogenesis and antiretroviral treatment ([4,14]). For a detailed mathematical analysis on such models, we refer to [15] and [5].

Fractional-order dynamics has recently been a focus of interest because of its appearance in physics, biology and engineering applications. There is a rich literature on the theoretical research of the fractional differential equations. The book of Podlubny [16] provides an overview to the basic theory of fractional differential equations. The monograph by Samko, Kilbas, and Marichev [17] contains remarkably comprehensive theory on the topic of fractional differential equations. Recently, much work has done on modeling the HIV infection with fractional-order derivatives [18-23] and [24]. A fractional-order model retains its memory, which is one of the main characteristic of the fractional-order derivative, while the features of immune response include memory. In [22] a fractional-order time-delay model is investigated which include three types of cells, namely healthy CD4+ T-cells, infected CD4+ T-cells and free HIV virus particles. In [19] the authors introduced fractional orders to the model of HIV-1 whose components are plasma densities of uninfected CD4+ T-cells, they use the generalized Euler method (GEM) and homotopy analysis method (HAM) to approximate the solution. Approximate solutions of fractional-order differential system for modeling human T-cell lymphotropic virus I (HTLV-I) infection of CD4+ T-cells is investigated in [25] using a multi-step generalized differential transform method. Bernstein operational matrices is applied to fractional order HIV model to approximate the solution [26]. Homotopy decomposition method is given in [27] to solve a system of fractional nonlinear differential equations that arise in the model for HIV infection of CD4+ T-cells. Recently, Huo etal. [28] studied a fractional-order HIV model to assess the impact of vaccines in a homosexual community. They have shown that the vaccinated reproduction below unity is not a threshold of HIV eradication when effectiveness and the dosage of the vaccines are low. A new critical threshold is derived in order to eradicate the HIV. Recently Pinto etal. [29] studied a fractional-order model for HIV infection that includes latently infected cells, macrophages and CTLs. In this paper we propose a finite difference implicit scheme to solve fractional-order HIV model containing four types of populations: uninfected CD4+ T-cells, latently infected, actively infected CD4+ T-cells and HIV virus particles. Further we analyze the dynamics of fractional-order HIV by investigating bifurcation points. The benefit of fractional-order systems is that they permit greater degrees of freedom in the model.

This paper is organized as follows: Section 2 includes generalized HIV model. In Section 3, we introduce implicit scheme for solving the fractional-order HIV-1 infection model. Dynamical behavior of generalized HIV model is investigated in Section 4. In Section 5, we present numerical simulations of the model and discuss the biological significance of the results. A conclusion is given in the last section.

2 Generalized HIV model

The Riemann-Liouville fractional integral Ia u of order a > 0 of u: R+ ^ R is defined by 1 f t

Ia u(t) = —1—- (t - s)a-1u(s) ds,

r(a) Jo

provided the expression on right hand side is defined. Here r denotes the Gamma function [30].

The Caputo fractional derivative Da u of order a of a continuous function u: R+ ^ R is defined by

i in(t- s)m-a-1u(m)(s) ds, m -1 < a < m,

Dau(t) = wmaJo( ) () -

m = a.

In particular, when 0 < a < 1, we have

r(1- a)

Dau(t)= w/ N f (t - s)-au (s) ds. ^ 1 "Jo

We generalize the integer-order model of target cells limited proposed by Perelson et al. [31] to fractional order 0 < a <1, which involves various types of cells: Let T, TL and TA denote the population of uninfected CD4+ T-cells, latently infected and actively infected CD4+ T-cells, respectively. The population of free virus particles is denoted by V. Interaction between these cells is given by the following system:

DaT = sa - iiaTT + ra T(1- T +TL+ta )-k?TV, t > 0,

1 T max 1

Da Tl = ka TV - ixaTTL - k2a Tl, t > 0, (1)

Da Ta = ka Tl - Ta , t > 0, ()

Da V = Niah Ta - kl VT - iaVV, t > 0.

Note that the units of the fractional differential equation are different, that is, fractional differential equations are expressed with respect to an intrinsic time variable that depends on a [32, 33] instead of the physical time. Thus a reformed parameter has been presented in model (1) to interpret the fractional order. Notice that when a ^ 1 the fractional HIV model (1) reduces to the classical HIV model.

The first two equations deal with the effects of HIV. Here sa, corresponds to the s, the source term, from the classical HIV model, ra corresponding to r represents the rate of growth for the CD4+ cell population and iaT is the analogon to the i T , the death rate of uninfected CD4+ T-cells. corresponding to k\ represents the rate at which CD4+ T-cell become infected by virus modeled by mass-action type of term. The dynamical behavior of actively infected T-cells is given in the third equation with ka corresponding to the k2 rate

at which latently infected cells convert to being actively infected and [ corresponding to [b represents the death rate per cell from the classical HIV model. The last equation models the free infectious virus population. Assume that an actively infected CD4+ T-cell produces N virus particles. The virus production rate is N times [. Free virus is expired at a rate k% VT by binding to uninfected CD4+ T-cells. Viral death from the body is represented by V.

3 Construction of implicit numerical scheme

In this section, we introduce an implicit numerical scheme using a finite difference approximation of the Caputo derivative. Numerous schemes have been developed for the numerical solution of the fractional differential equations [34]. A class of the fractional multi-step method is proposed by Lubich [35] and Galeone and Garrappa [36], the fractional Adams method is proposed by Diethelm et al. [37] and Odibat and Momani [3] and Grunwald-Letnikov approximation based on the Grunwald-Letnikov definition of the fractional derivative is addressed in [38] and discussed the analysis of convergence and stability. Baker [39] introduced some numerical methods for the Volterra integral and integro-differential equations.

In this paper, we implement the Ll-scheme to approximate the Caputo fractional derivative, which was independently developed and analyzed in [40] and [41]. The Ll-scheme is based on a piecewise linear approximation to the fractional derivative. We are in favor of the Ll-scheme, because this scheme is derived easily and the coefficients of this scheme have good properties e.g. the representation of the coefficients is simple. The Ll-scheme has been extensively used in practice and currently it is one of the most efficient numerical methods for solving the time fractional differential equations due to its ease of implementation.

To discretize Daf (t) based on the Ll-scheme, first we defined the temporal size t and tn means m. Therefore we discretize the Caputo derivative by a finite difference method,

Daf (tn) =

r(l- a)

f (tn - ,)-f) ds

l n-l r t+l

V r\t„ - s)f+lf) ds + R

r(l-a)^ J,

\ ) j=o Jtj

r(2 - a)

Y,bJf (tj)+Rn,

-(nl-a -(n - l)l-a), j = 0,

(n - j + l)l-a - 2(n - j)l-a + (n - j - l)l-a, j = l, 2,..., n -1,

l, j = n.

V^ (tn) =

r(2 - a)

Y,bjf (tj), n = l, 2,..., N,

Algorithm 1 A modified Jacobi matrix on Newton iteration method

1: Given Witer, e, let Nd = 1 2: if k > Niter and ||f(xM)\\ > ||f(xk)||

3: Nd :=2 Nd 4: end if

5: if det(J(xk)) < e, %e is a sufficiently small positive number

6: J(xk) := Nd ■ (J(xk) +1), %l is an identity matrix 7: else

8: J(xk):= Nd ■ J(xk) 9: end if

then by [41], we have, for n = 1,2,..., N,

\Rnt\ = |£>V (in) - Vf (in) | < C max

d2f (t) di2

for some C > 0. The implicit numerical scheme for fractional-order system (1) with finite difference approximation of the Caputo derivative is

i Tt+Tt+Tt

Va Tn = sa - /aTTn + ra Tn(1 - 1 +TL+±A) - № Tn Vn,

L 1 1 ma^ 1

VaTl = k1 Tn Vn - fiTTl - Tl, va TA = ka ti - /4 TA, VaVn = N/l Tn - kf VnTn - /iaVVn,

that is,

1 Tn + Tn + Tn

Tn + j bn T> - Ta r(2 - a)(sl - /aTTn + ra Tn(1 - 1 +t1Jl+1a ) - k? TnVn) = 0, TL + E^-)1 b-Tl - Tlr(2 - a)kTnVn - /H - TnL) = 0,

TA + EI-o1 % T^ - Tl r (2 - a)(ka Tin - / TA) = 0,

Vn + E^-,1 b"V> - Tlr(2 - a)(N/ TAn - k? VnTn - /aVVn) = 0.

In order to solve the above nonlinear equations, we choose the Newton iteration method

xk+1 = xk -/-1 (xk )f (xk), k = 0,1,2,..., (6)

the initial values are given by x0 = (T0, T), T), V0). Furthermore, to ensure the convergence of the Newton iteration method and avoid the Jacobian matrix /(xk) to be nonsingular, we improve the above method by the Algorithm 1.

4 Dynamical behavior of generalized HIV model

In this section we will analyze the dynamics of the generalized HIV model and examine the effect of the fractional order on the HIV dynamics. For this we find the equilibrium points of (1) by solving the following system:

' sa - /TT + ra T(1 - T+1L+TA )-k? TV = 0,

T Tmax 1

ka TV - /TTL - ka Tl = 0, kk? Tl - =0, N/I\TA - klVT - /VV = 0.

We find that system (1) has two equilibrium points: the uninfected equilibrium points Ef = (T0, 0,0,0) and the (positive) infected equilibrium points Ea = (T, TL, TA, V) where

^/p2 + 4sa

¡¿Vh

ka (Nka - k3y

¡vv ka/xavV

tl = TTTa—t , ta =

Nka - k^ A ¡a (Nka - k3)'

p = ra - ¡T, and Y = /Tmax,

V = sa k|+ pk4 ¡V - Y (¡V)2

ka^av (k4 + k5^V ) ,

k4 = kf

ka+¡4

-1, k3 = ka + ¡T, k5 = 1 +

Theorem 4.1 Tfe uninfected equilibrium point Ea = (To, 0,0,0) is asymptotically stable if and only if

¡V > ¡crit =

kg kg ToN - kg k3To

Proof The Jacobian matrix Ji evaluated at E0a for the system (1) is given by

-a -Y T0 -Y T0 -ka ta

0 -k3 0 ka T0

0 ka -¡b 0

\ 0 0 N ¡iab -ke

where ke = kaT0 + ¡¡V, a = -p + 2T0y. The associated characteristic equation is given by

(X + a)[ (X + ¡a) (X + k3)(X + k6) - ka ka TNtf,] = 0. Hence

X1 = -a < 0,

since substituting the value of T0, we get a = ^p2 + 4sy > 0. Dividing the characteristic equation by (X + a), we obtain

X3+ A1X2 + A2X + A3 = 0,

Ai = ¡a + k3 + ke > 0, A 2 = ka ke + ¡a (k3 + ke) > 0, A3 = ¡1 (k3ke- kia k2 T0N).

Using the definition of we have

A3 = rf {hk6- ka T0N)

= ^ (k3(k? t0+^ - ka ka t0N)

'ki ka ToN - ka k3To

= l4 hy¿V -

= ¿1 kf{¿V - ¿crit).

Eigenvalues have a negative real part if and only if the following conditions of the Routh-Hurwitz criteria are satisfied:

Ai, A3>0 and AiA2- A3 > 0.

Here Ai is positive, being a sum of positive terms,

A1A2 - A3 = ¿a (k3 + k6) + ¿a (kf + k2 + 2k3k6 + k?k? ToN) + k3k6(k3 + k6) > 0.

A3 > 0 if and only if

lxV > !xcrit'

Remark 4.1 ¿V is a bifurcation parameter, transcritical bifurcation occurs as ¿V passes the point ¿rit (see Figure i). When

> rfcrit =

kf kg ToN - kf k3To

the uninfected equilibrium point Eo is stable and the infected equilibrium point Ef does not exist (this case is unphysical). When ¡xV = rfOrit< the uninfected equilibrium point is neutrally stable. When ¡xV < Eo becomes unstable.

Remark 4.2 Following a similar analysis, we see that a transcritical bifurcation occurs at Nrit, where

k3(rfV + kf To)

crit kl k% To '

If N < N?rit; the uninfected equilibrium point E^ is stable and the infected equilibrium point E^ does not exist (this case is unphysical). The uninfected equilibrium point is neutrally stable at N = N?rit. Eo becomes unstable when N > Ncarit.

The Jacobian matrix /2 evaluated at Ef for the system (1) is given by

i -a -y T -yT -kf T\ 0

-kf T -k3

\-kf T o N-k6 /

where a = -p + 2y T + y T + k1 V. The characteristic equation associated to the Jacobian matrix is given by

X4 + B3k3 + B2X2 + Bik + B0 = 0, (8)

with Bo = kaVT[kZvab(Nk% - ka) + y^v(ka + ¡%)], Bi = a[k3k6 + ^(k3 + ke)j + KVT x [y(¡V + ka + Iiab) - kia(k3 + ¡¡¡)], B2 = a(k3 + k6 + ¡4) + ^g(k3 + k6) + k3k4 + ka TV(y - kia), B3 = a + k3 + k6 + ¡a.

By the Routh-Hurwitz criteria for the stability of fractional-order systems [42], we have the following result for the stability of infected equilibrium.

Theorem 4.2 The infected equilibrium point E^ = (T, TL, TA, V) is locally asymptotically stable if the coefficients of the characteristic polynomial (8) evaluated atEf satisfy

Bo > 0, Bi>0, B3 > 0 and B1B2B3 > B2 + B2Bo, (9)

for all a e (0,1].

5 Numerical simulations

In this section, numerical simulations are provided to verify the theoretical results established in the previous sections. In the following we will monitor the effect of varying viral death rate and varying number of new virions produced by infected CD4+ T-cells on the dynamical behavior of the fractional model, respectively. The parameters are chosen as in Table 1, unless otherwise stated.

Table 1 HIV model parameters

¿a Ib I'V kg kg ra N Tmax s"

0.02a 0.24a 2.4a (2.4 x 10-5)a (3 x 10-3)a 0.03 Varies 1,500 10a

IV (Days)

Figure 1 Bifurcation diagram (viral death rate), N = 1,000.

1000 900

1800 I-

4 700 o

| 600 ®

Jj 920 3

— O=0.99

-—№=0.97

- - a=0.95

0 100 200 300 400 500 600 700 BOO 900 1000 t-time (Days)

0 100 200 300 400 500 600 700 800 900 1000 t-time (Days)

ct=0.99 a=0.97 a=0.95

0 100 200 300 400 500 600 700 800 900 1000 t-time (Days)

0 100 200 300 400 500 600 700 800 900 1000 t-time (Days)

- a=0.

- a=0.97

- - 0=0.95

"t----1-----1----1

0 100 200 300 400 500 600 700 800 900 1000 1-time (Days)

0 100 200 300 400 500 600 700 800 900 1000 t-time (Days)

1300 >

£200 100

0 100 200 300 400 500 600 700 800 900 1000 t-time (Days)

— a=0.99 - a=0.97 a=0.95

0 100 200 300 400 500 600 700 800 900 1000 t-time (Days)

Figure 2 Solution of fractional HIV model with fl^ = 2.4a (left panel) and fl^ = 7.4a (right panel) for

Uninfected equilibria for different fractional orders are as follows:

' E? = (978,0,0,0), ¿arit = 3.4364 for a = 0.99, E? = (936,0,0,0), ¿aait = 4.2059 for a = 0.97, E? = (986,0,0,0), ¿aait = 5.1488 for a = 0.95.

By Remark 4.1, Ea is asymptotically stable when ¿V = 7.4a > ¿arit and Ea become unstable as ¿V = 2.4a is reduced. Transcritical bifurcation occurs at the point ¿arit see Figure i.

500 1000

t-time (Days)

500 1000

1-tims (Days)

120 r = 100

a=0.99 o=0.97 o=0.95

100 200 300 400 500 600 700 800 900 1000 t-tlms (Days)

ci=Q.99 ct=0.97 a=0.95

100 200 300 400 500 600 700 800 900 1000 t-time (Days)

a> o £

— a=0 99

- -ct=0,97

100 200 300 400 500 600 700 800 900 1000 t-time (Days)

— a=0.99

— -o=0.97

— - o=0.95

100 200 300 400 500 600 700 800 900 1000 t-time (Days)

o=0.97 - o=Q ,95

— o=0.99

- o=0.97 -- o=0.95

100 200 300 400 500 600 700 800 900 1000

t-time (Days)

100 200 300 400 500 600 700 800 900 1000

t-time (Days)

Figure 3 Solution of fractional HIV model with N = 600 (left panel) and N = 500 (right panel) for a =0.99,0.97,0.95 with initial condition (7°, 7"°, 7°, V0) = (1,000,0.01,0.1,0.0001).

The approximate solution for a = 0.99,0.97,0.95 are displayed in Figure 2, which are in good agreement with the analytical result.

Next we will simulate the fractional system by varying the values of the parameter N and choosing other values from Table 1. By calculation we can obtain

E^ = (978,0,0,0), Nrit = 694 for a = 0.99, E^ = (936,0,0,0), Ncarit = 559 for a = 0.97, Ea = (986,0,0,0), NCTit = 450 for a = 0.95.

According to Remark 4.2, E°-99 is stable but E°-97 and E°-95 are unstable when N = 600. Next reducing N = 500 gives N < Nc°rft9 and N < Nc°n9t7 so E0'99 and E0'97 become stable; however, E°'95 is unstable as shown in Figure 3.

For a = 0.97 the coefficients of characteristic polynomial evaluated at infected equilibrium EH = (520,236,3,359) are

B0 = 2.6259 x 10-4 >0, B1 = 0.0197 >0, B3 = 2.6610 > 0, B1B2B3 = 0.0353 > B2 + B33B0 = 0.0022.

Hence by Theorem 4.2, E°-97 is locally asymptotically stable as shown in Figure 4.

The decline in the number of CD4+ T-cells in peripheral blood is used in a clinical setting as indicators of the disease stage. Decreasing fractional order gives rise to larger amounts of CD4+ T-cell depletion as shown in Figure 5.

Figure 5 Effect of fractional derivative on T-cell depletion for N = 900,1,200 other parameters are given in Table 1 with initial condition (1,000,0,0,0.001).

6 Conclusion

In this paper, we have investigated a fractional-order HIV model, as a generalization of an integer-order model. The advantage of the generalized model is that the fractional-order system possesses memory, which belongs to the main features of the immune response. An implicit numerical scheme has been proposed for the fractional-order HIV model using a finite difference approximation of the Caputo derivative. We showed that the system possesses two equilibrium points and derived the analytical condition for the stability of uninfected and infected equilibrium points. Further we analyzed the influence of the fractional derivative on the dynamics of system. In many AIDS patients the T-cell level can be depleted <200 mm-3 level, on the other hand an integer-order model cannot model this fact using the parameter values in Table I. However, with the additional degree of freedom of the fractional derivative we can obtain the low CD4+ cell counts seen during the late stages of the disease; see Figure 5. We established that ¿V and N are bifurcation parameters, transcritical bifurcation occurs as ¿V and N passes the point ¿rit and N^, respectively. We have shown that the disease can be eradicated by increasing the viral death rate greater than ¿nt as shown in Figure 2. Another parameter that could play a vital role to control the HIV virus is the number of new virions produced by infected CD4+ T-cells. We found that the virus can be eliminated if N is less than Narit, that is, HIV infection will not be sustained if infected cells die without producing an adequate number of viral progeny as demonstrated in Figure 3. Mathematical models can help physicians to choose an optimal dosage and check the effects of their therapeutic action. HIV progression has many variations from patient to patient, which is difficult to capture by an integer-order derivative. The fractional derivative can be varied to best fit the real data according to the progression of different HIV patients. Thus a more reliable model can be obtained by choosing the relevant fractional index according to real data. Consequently, the clinician can recommend administration of drugs or treatment to each individual patient by using the information from the generalized model with the most relevant fractional index.

Competing interests

The authors declare that they have no competing interests. Authors' contributions

All authors contributed equally to this work. They all read and approved the final version of the manuscript. Author details

1LSEC, ICMSEC, Academy of Mathematics and Systems Science, Chinese Academy of Sciences, Beijing, 100190, China. 2Department of Mathematics, Cankaya University, Ankara, 06530, Turkey. 3School of Mathematics and Computational Science, Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Xiangtan University, Xiangtan, 411105, China. institute of Space Sciences, Magurele-Bucharest, 077125, Romania. 5School of Mathematical Sciences, University of Chinese Academy of Sciences, Beijing, 100049, China. 6COMSATS Institute of Information Technology, Lahore, Pakistan.

Acknowledgements

The authors would like to thank the anonymous referees for their valuable suggestions. This research is supported by the State Key Program of National Natural Science Foundation of China (Grant No. 11631013), the National Natural Science Foundation of China (Grant No. 11371357) and by the ITER-China Program (Grant No. 2014GB124005). The third author extend his appreciation to the National Natural Science Foundation of China (No. 11601460) and the Research Foundation of Education Commission of Hunan Province of China (No. 16C1540) for supporting this research work.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Received: 30 December 2016 Accepted: 17 March 2017 Published online: 29 March 2017

References

1. Nowak, MA, Bangham, CR: Population dynamics of immune responses to persistent viruses. Science 272, 74-79 (1996)

2. Nowak, MA, Bonhoeffer, S, Shaw, GM, May, RM: Anti-viraldrug treatment: dynamics of resistance in free virus and infected cell populations. J. Theor. Biol. 184, 203-217 (1997)

3. Odibat, Z, Momani, S: An algorithm for the numericalsolution of differentialequations of fractionalorder. J. Appl. Math. Inform. 26,15-27 (2008)

4. Perelson, AS: Modelling viraland immune system dynamics. Nat. Rev. Immunol. 2, 28-36 (2002)

5. Perelson, AS, Nelson, PW: Mathematicalanalysis of HIV-1 dynamics in vivo. SIAM Rev. 41, 3-44 (1999)

6. Perelson, AS, Essunger, P, Cao, Y, Vesanen, M, Hurley, A, Saksela, K, Markowitz, M, Ho, DD: Decay characteristics of HIV-1-infected compartments during combination therapy. Nature 387,188-191 (1997)

7. Perelson, AS, Essunger, P, Ho, DD: Dynamics of HIV-1 and CD4+ lymphocytes in vivo. AIDS 11, S17-S24 (1997)

8. Perelson, AS, Neumann, AU, Markowitz, M, Leonard, JM, Ho, DD: HIV-1 dynamics in vivo: virion clearance rate, infected celllife-span, and viralgeneration time. Science 271,1582-1586(1996)

9. Silva, CJ, Torres, DFM: ATB-HIV/AIDS coinfection modeland optimalcontroltreatment. Discrete Contin. Dyn. Syst. 35(9), 4639-4663 (2015)

10. Rocha, D, Silva, CJ, Torres, DFM: Stability and optimalcontrolofa delayed HIV model. Math. Methods Appl. Sci. (2016). doi:10.1002/mma.4207

11. Luo, J, Wang, W, Chen, H, Fu, R: Bifurcations of a mathematicalmodelfor HIV dynamics. J. Math. Anal. Appl. 434, 837-857(2016)

12. Wang, Y, Liu, J, Liu, L: Viraldynamics of an HIV modelwith latent infection incorporating antiretroviraltherapy. Adv. Differ. Equ. 2016, Article ID 225 (2016)

13. Rong, L, Perelson, AS: Modeling HIV persistence, the latent reservoir, and viralblips. J. Theor. Biol. 260, 308-331 (2009)

14. Finzi, D, Siliciano, R: Viraldynamics in HIV-1 infection. Cell 93,665-671 (1998)

15. Kirschner, DE: Using mathematics to understand HIV immune dynamics. Not. Am. Math. Soc. 43, 191-202 (1996)

16. Podlubny, I: Fractional DifferentialEuations. Acdemic Press, San Diego (1999)

17. Samko, SG, Kilbas, AA, Marichev, OI: Fractional Integrals and Derivatives: Theory and Applications. Gordon and Breach Science Publishers, Yverdon (1993)

18. Arafa, AAM, Rida, SZ, Khalil, M: Fractionalmodeling dynamics of HIV and CD4+ T-cells during primary infection. Nonlinear Biomed. Phys. 6, Article ID 1 (2012)

19. Arafa, AAM, Rida, SZ, Khalil, M: The effect of anti-viraldrug treatment of human immunodeficiency virus type 1 (HIV-1) described by a fractionalorder model. Appl. Math. Model. 37, 2189-2196 (2013)

20. Gökdogana, A, Yildiri m, A, Merdana, M: Solving a fractionalorder modelof HIV infection of CD4+ T cells. Math. Comput. Model. 54, 2132-2138 (2011)

21. Kou, CH, Yan, Y, Liu, J: Stability analysis for fractionaldifferentialequations and their applications in the models of HIV-1 infection. Comput. Model. Eng. Sci. 39, 301-317 (2009)

22. Liu, Z, Lu, P: Stability analysis for HIV infection of CD4+ T-cells by a fractional differentialtime-delay modelwith cure rate. Adv. Differ. Equ. 2014, Article ID 298 (2014)

23. Pinto, CMA, Carvalho, ARM: Fractionalmodeling of typicalstages in HIV epidemics with drug-resistance. Prog. Fract. Differ. Appl. 1(2), 1 1 1-122 (2015)

24. Yan, Y, Kou, C: Stability analysis for a fractionaldifferentialmodelof HIV infection of CD4+ T-cells with time delay. Math. Comput. Simul. 82, 1572-1585 (2012)

25. Ertürk, VS, Odibat, ZM, Momanic, S: An approximate solution of a fractionalorder differentialequation modelof human T-celllymphotropic virus I (HTLV-I) infection of CD4+T-cells. Comput. Math. Appl. 62, 996-1002(2011)

26. Alipour, M, Arshad, S, Baleanu, D: Numericaland bifurcations analysis for multi-order fractionalmodelof HIV infection of CD4+T-cells. Sci. Bull. "Politeh." Univ. Buchar., Ser. A, Appl. Math. Phys. 78(4), 243-258 (2016)

27. Atangana, A, Alabaraoye, E: Solving a system offractional partialdifferentialequations arising in the modelof HIV infection of CD4+ cells and attractor one-dimensionalKeller-Segelequations. Adv. Differ. Equ. 2013, Article ID 94 (2013)

28. Huo, J, Zhao, H, Zhu, L: The effect of vaccines on backward bifurcation in a fractionalorder HIV model. Nonlinear Anal., RealWorld Appl. 26, 289-305 (2015)

29. Pinto, CMA, Carvalho, ARM: A latency fractionalorder modelfor HIV dynamics. J. Comput. Appl. Math. 312, 240-256 (2017)

30. Abramowitz, M, Stegun, IA (ed.): Handbook of MathematicalFunctions: With Formulas, Graphs, and Mathematical Tables. NationalBureau of Standards Applied Mathematics Series, vol. 55. U.S. Government Printing Office, Washington (1964)

31. Perelson, AS, Kirschner, DE, De Boer, R: Dynamics of HIV infection of CD4 T-cells. Math. Biosci. 114, 81-125 (1993)

32. Arenas, AJ, Gonzalez-Parrab, G, Chen-Charpentierc, BM: Construction of nonstandard finite difference schemes for the SI and SIR epidemic models of fractionalorder. Math. Comput. Simul. 121,48-63 (2016)

33. Diethelm, K: A fractionalcalculus based modelfor the simulation of an outbreak of dengue fever. Nonlinear Dyn. 71, 613-619(2013)

34. Weilbeer, M: Efficient numericalmethods for fractionaldifferentialequations and their analyticalbackground. Ph.D. thesis,Technischen Universität Braunschweig (2005)

35. Lubich, C: Discretized fractionalcalculus. SIAM J. Math. Anal. 17, 704-719 (1986)

36. Galeone, L, Garrappa, R: On multistep methods for differentialequations of fractionalorder. Mediterr. J. Math. 3, 565-580 (2006)

37. Diethelm, K, Ford, NJ, Freed, AD: A predictor-corrector approach for the numericalsolution of fractionaldifferential equations. Nonlinear Dyn. 29, 3-22 (2002)

38. Scherer, R, Kalla, SL, Tang, Y, Huang, J: The Grünwald-Letnikov method for fractionaldifferentialequations. Comput. Math. Appl. 62, 902-917(2011)

39. Baker, CTH: A perspective on the numericaltreatment of Volterra equations. J. Comput. Appl. Math. 125, 217-249 (2000)

40. Lin, Y, Xu, C: Finite difference/spectralapproximations for the time-fractionaldiffusion equation. J. Comput. Phys. 225(2), 1533-1552 (2007)

41. Sun, Z, Wu, X: A fully discrete difference scheme for a diffusion-wave system. Appl. Numer. Math. 56,193-209 (2006)

42. Ahmed, E, El-Sayed, AMA, El-Saka, HAA: On some Routh-Hurwitz conditions for fractionalorder differentialequations and their applications in Lorenz, Rossler, Chua and Chen systems. Phys. Lett. A 358,1 -4 (2006)

Submit your manuscript to a SpringerOpen journal and benefit from:

► Convenient online submission

► Rigorous peer review

► Immediate publication on acceptance

► Open access: articles freely available online

► High visibility within the field

► Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com