Journal of the Egyptian Mathematical Society (2013) xxx, xxx-xxx
Egyptian Mathematical Society Journal of the Egyptian Mathematical Society
www.etms-eg.org www.elsevier.com/locate/joems
A fractional-order model of HIV infection with drug therapy effect
A.A.M. Arafa a *, S.Z. Rida b, M. Khalilc
a Department of Mathematics & Computer Science, Faculty of Science, Port Said University, Port Said, Egypt b Department of Mathematics, Faculty of Science, South Valley University, Qena, Egypt c Department of Mathematics, Faculty of Engineering, October University For Modern Science and Arts, MSA University, 6th Oct. City, Giza, Egypt
Received 25 April 2013; revised 5 November 2013; accepted 11 November 2013
KEYWORDS
Fractional calculus; Euler method; HIV model
Abstract In this paper, the fractional-order model that describes HIV infection of CD4+ T cells with therapy effect is given. Generalized Euler Method (GEM) is employed to get numerical solution of such problem. The fractional derivatives are described in the Caputo sense.
2000 MATHEMATICS SUBJECT CLASSIFICATION: 34A08; 92Bxx
© 2013 Production and hosting by Elsevier B.V. on behalf of Egyptian Mathematical Society.
1. Introduction
HIV is a retrovirus that targets the CD4+ T lymphocytes, which are the most abundant white blood cells of the immune system [2,3,26,27]. Although HIV infects other cells also, it wreaks the most havoc on the CD4+ T cells by causing their decline and destruction, thus decreasing the resistance of the immune system [20-25]. Mathematical models have been proven valuable in understanding the dynamics of HIV infection [5,6,15,18,19]. Anti-viral drug treatment for human immunodeficiency virus (HIV) infections causes rapid reduction in
* Corresponding author.
E-mail addresses: anaszi2@yahoo.com (A.A.M. Arafa), szagloul@
yahoo.com (S.Z. Rida), m_kh1512@yahoo.com (M. Khalil).
Peer review under responsibility of Egyptian Mathematical Society.
plasma virus load. Reverse transcriptase inhibitors (RTIs) are a class of antiretroviral drug used to treat HIV infection [23,24]. RTIs inhibit [23] activity of reverse transcriptase, a viral DNA polymerase enzyme that retroviruses need to reproduce. In HIV-1, reverse transcriptase inhibitors block the infection of uninfected cells [25,26]. The authors in [25], assumed that The drug may not be 100% effective, so a part of infected cells in pre-RT class will revert back to uninfected class and the remaining will progress to complete reverse transcription and become productively infected and then produce virus. They developed a mathematical model for primary infection with RT inhibitor under the above mentioned assumptions. In view of this, the following model is a generalization of the model presented in [25]:
Da(T) — s - kLT- iT + (ge + b)I Da(I) — kLT -(i + e + b)I Da(V) — (1 - g)eI- dV D"(L) — NSV- cL where 0 < a 6 1.
1110-256X © 2013 Production and hosting by Elsevier B.V. on behalf of Egyptian Mathematical Society. http://dx.doi.Org/10.1016/j.joems.2013.11.001
In the above model, T represents density of susceptible CD4+ T cells, I represents density of infected CD4+ T cells before reverse transcription (pre-RT class), V represents density of infected CD4+ T cells in which reverse transcription is completed (post-RT class) and they are capable of producing virus. V is virus density. After infection, infected cells progress to pre-RT after infection, infected cells progress to pre-RT class I and then they leave pre-RT class at a rate e to productively infected (post-RT) class. These infected cells are capable of producing virus particles. In view of the above discussion, we consider that due to the presence of RT inhibitor a fraction of cells g el in pre-RT class reverts back to uninfected class and remaining (1 — g)eI proceeds to post-RT class and become productively infected, where 0 < g < 1 is the efficacy of RT inhibitor. The parameter s is the inflows rate of CD4+ T cells, l its natural death rate, k is interaction infection rate of CD4+ T cells, h is death rate of infected cells, e is the transition rate from pre-RT infected CD4+ T cells class to productively infected class (post-RT), b is the reverting rate of infected cells to uninfected class due to non-completion of reverse transcription [23,25]. The parameter d represents death rate of actively infected cells and includes the possibility of death by bursting of infected T cells, c is the clearance rate of virus and N is the total number of viral particles produced by an infected cell. The rest of the paper is organized as follows. In Section 2, the idea of the fractional calculus is presented, while Generalized Euler method (GEM) is presented in Section 3.Non-nega-tive solutions of the presented system are discussed in Section
4.The equilibrium points and stability are discussed in Section
5. Section 6 is devoted for the numerical results of the presented problem (1).
terpart [1]. Hence, we propose a system of FODE for modeling HIV infection. We first give the definition of fractional-order integration and fractional-order differentiation [10]. For the concept of fractional derivative, we will adopt Caputo's definition, which is a modification of the Riemann-Liouville definition and has the advantage of dealing properly with initial value problems.
Definition 2.1. Riemann-Liouville fractional integration of order a is defined as:
fx) (x — t)a—1Mdt, a > 0, x > 0
C (a) Jo
J0f(x) =Ax)
Definition 2.2. Riemann-Liouville and Caputo fractional derivatives of order a can be defined respectively as follows:
Daf(x) = Dm(Jm—af(x)), Dfx) = Jm—a(Dmf(x)) where
m — 1 < a 6 m, m 2 N
(2) (3)
Properties of the operator J can be found in [4,5], we mention only the following:
(1) JaJbf(x) = Ja+bf(x)
(2) JaJbf(x) = JbJaf(x)
(3) Jatc =
rêrTTÏÏf+' a>c> —1, '>0 (6)
2. Fractional calculus
3. Generalized Euler Method (GEM)
Fractional calculus (FC) has been extensively applied in many fields [4,7]. Many mathematicians and applied researchers have tried to model real processes using the fractional calculus. Jesus, Machado and Cunha analyzed the fractional-order dynamics in botanical electrical impedances [11,12]. Petrovic, Spasic and Atanackovic developed a fractional-order mathematical model of a human root dentin [20]. In biology, it has been deduced that the membranes of cells of biological organism have fractional-order electrical conductance [4] and then are classified in groups of non-integer order models. Fractional derivatives embody essential features of cell rheological behavior and have enjoyed greatest success in the field of rheology [7]. Also, it has been shown that mathematical modeling by fractional ordinary differential equations (FODE) has more advantages than classical integer order modeling [9]. The major reason of using is that fractional differential equations are naturally related to systems with memory which exists in most biological and systems [1-4]. In other words, Calculating time-fractional derivative of a function f(t) at some time t = tx requires all the previous history, i.e. all f(t) from t = 0 to t = t1. Also, using fractional-order differential equations can help us to reduce the errors arising from the neglected parameters in modeling real life phenomena [6]. FODE are naturally related to systems with memory which exists in most biological systems. Also, they are closely related to fractals, which are abundant in biological systems [8]. Fractional-order differential equations are, at least, as stable as their integer order coun-
Most non-linear fractional differential equations do not have analytic solutions [16,17], so approximations and numerical techniques must be used. The decomposition method (ADM) and the variational iteration method (VIM) are relatively new approaches to provide an analytical approximate solution to linear and non-linear problems, and they are particularly valuable as tools for scientists and applied mathematicians, because they provide immediate and visible symbolic terms of analytic solutions, as well as numerical approximate solutions to both linear and non-linear differential equations. In recent years, the application of the ADM, VIM, in linear and nonlinear problems has been developed. On the other hand, these methods are effective for small time, i.e. t 1 [9,10], however the standard homotopy perturbation method (HPM) [14] cannot solve the problem for larger time and in fact the solution of the chaotic system using HPM is an open problem. Nevertheless by chance, there are cases at which these methods give good approximation for a large range of time (t). A few numerical methods for fractional differential equations have been presented in the literature. However many of these methods are used for very specific types of differential equations, often just linear equations or even smaller classes. In [16,17], Odibat and Momani derived the generalized Euler's method that we have developed for the numerical solution of initial value problems with Caputo derivatives. The method is a generalization of the classical Euler's method. Consider the initial value problem
Dly(t) —f(t;y(t)), y(0) —y0; 0 < a 6 1; t > 0
Let [0, a] be the interval over which we want to find the solution of the problem (7). In actuality, we will not find a function y(t) that satisfies the initial value problem (7). Instead, a set of points {(tj, y(tj)} is generated, and the points are used for our approximation. For convenience we subdivide the interval [0, a] into k subintervals [tj, tj +1] of equal width h = a/k by using the nodes tj = jh, for j = 0, 1, ..., k. Assume that y(t),Daty(t), and D2tay(t) are continuous on [0, a] and use the generalized Taylor's formula to expand y(t) about t = t0 = 0. For each value t there is a value c1 so that
y(t)—y(t0) + {D"-y(t))(t0 )
C(a + 1)
+ Dymc!)
C(2a + 1)
When (Daty(t))(t0)= f(t0,y(t0)) and h = tx are substituted into Eq. (8), the result is an expression for y(t1):
y(h)— y(h)+f(t0, y(t0))
ha h h + D^ma) h
C(a + 1)
C(2a + 1)
If the step size h is chosen small enough, then we may neglect the second-order term (involving h2a) and get
y(h)—y(h)-
C(a + 1)
f(t0; y(t0))
The process is repeated and generates a sequence of points that approximates the solution y(t). The general
formula
y(tj+\) — y(tj)
for generalized is
Euler's method (GEM) when
tj + h
ftj ; y(tj))
C(a + 1)
for j = 0, 1, ..., k — 1. It is clear that if a = 1, then the generalized Euler's method (9) reduces to the classical Euler's method.
4. Non-negative solutions
Denote = {X 2 R4 : X P 0g, and X(t) = (T, I, V, L)T. For the proof of the non-negative solution, consider the following theorem and corollary:
Theorem 1 (Generalized mean value theorem). Let f(x) 2 C(0, a] and Daf(x) 2 C(0, a], for 0 < a 6 1. Then we have
f(x) ~f(0+) + (Daf)(n)(x)a C(a)
with 0 6 n 6 x, "x 2 (0, a]. Proof. Proof is given in [16]. □
Corollary 1.1. Suppose that f(x) 2 C[0, a] and Daf(x) 2 C(0, a] for 0 < a 6 1. It is clear from theorem 1.1 that if Daf(x) P 0, "x 2 (0, a), then f(x) is non-decreasing and if Daf(x) 6 0, "x 2 (0, a) then f(x) is non-increasing for all x 2 [0, a].
Proof. This is clear from Theorem 1.1 [6]. □
Theorem 2. There is a unique solution X(t) = (T,I, V, L)Tfor (1) at t P 0 and the solution will remain in R4+.
Proof. The existence and uniqueness of the solution of the initial value problem (2.2) in (0, 1) can be obtained from [13, Theorem 3.1 and Remark 3.2]. Now we will show that R* is positively invariant domain:
DaT\T=0 = S + (ge + b)I P 0 DaI\I=0 = kLT P 0 Da V\ v=0 = (1 — g) eI P 0 DaL\L=0 = NdV P 0
From Corollary 1.1, the solution will remain in R4 . □
Ei = , 0,0, 0 ) and E2 = (T, I, V, L)
where T= (1i +e
NKe(i — g)
s — Tß
V= (I — g)eTl
-_NT*è
A sufficient condition for the local asymptotic stability of the equilibrium points is that the eigenvalues kt of the Jacobian matrix of Ei and the Jacobian matrix of E2 satisfy the condition | arg | > a This confirms that fractional-order differential equations are, at least, as stable as their integer order counterpart.
6. Numerical results
5. Equilibrium points and stability
The authors in [25] deduced the equilibrium Points of the integer system of the given system (1), i.e. when a = 1. To evaluate the equilibrium points of the fractional-order system (1), let
Da(T) = 0 Da(I) = 0 Da(V) = 0 Da(L) = 0
Then we will have the same results in [20] as follows
We will solve the system (i) by using (GEM). Consider that s = i0, k = 0.000024, d = 0.26, and c = 2.4. We choose l = 0.01 and ii = 0.015 (since death rate of cells with viral particle will be slightly higher than those of uninfected cells) and N = 1000. The drug efficacy g, e and b vary with initial condition T(0) = 300, 1(0) = 10, V(0) = 10, L(0) = 10.
7. Conclusion
In this paper we employed the Generalized Euler method (GEM) as a reasonable basis for studying the solution of human T-cell lymphotropic virus (HIV-I) infection of
Fig. 2 The densities of the susceptible CD4+ T cells T(t), infected CD4+ T cells I(t) in pre-RT class, density of infected CD4+ T cells V(t) in post-RT class, and the virus density L(t) when a = i. The gray line (e = 0.5), the dashed line (e = 0.4), black solid line (e = 0.3).
Fig. 3 The densities of the susceptible CD4+ T cells T(t), infected CD4+ T cells I(t) in pre-RT class, density of infected CD4+ T cells V(t) in post-RT class, and the virus density L(t) when a = 1. The gray line (b = 0), the dashed line (b = 0.05), black solid line (b = 0.075), the dotted line (b = 0.1), keeping e = 0.4 and g = 0.6.
Fig. 4 The densities of the susceptible CD4+ T cells T(t), infected CD4+ T cells I(t) in pre-RT class, density of infected CD4+ T cells V(t) in post-RT class, and the virus density L(t) for a =1 (the gray line) a = 0.99, (the dashed line), a = 0.98 (the black solid line).
CD4+ T cells. As g increases from 0.6 to 0.8 the viral level decreases and in case when g = 0.9 it approaches to zero (see Fig. 1), while the level of CD4+ T cells increases with increase in g. variation in e does not have much difference in viral level, though increase in a increases the viral level (see Fig. 2). For b = 0 to 0.1 again for same set of parameters as above. Keeping g = 0.6 and e = 0.4. The results are plotted in Fig. 3. It can be easily seen that increase in b shows increase in CD4+ T cells and decrease in viral level. The concentration of susceptible CD4+ T cells T(t), infected CD4+ T cells I(t), and free HIV virus particles V(t) in the blood have been obtained, therefore when a fi 1 the solution of the fractional model (1) Ta(t), Ia(t), Va(t), La(t), reduce to the standard solution T(t), I(t), V(t), L(t), (see Fig. 4). Finally, the recent appearance of fractional differential equations as models in some fields of applied mathematics makes it necessary to investigate methods of solution for such equations(analytical and numerical) and we hope that this work is a step in this direction.
Acknowledgment
The author would like to thank the anonymous reviewers very much for their valuable suggestions on improving this paper.
References
[1] E. Ahmed, A.M.A. El-Sayed, H.A.A. El-Saka, Equilibrium points, stability and numerical solutions of fractional-order predator-prey and rabies models, J. Math. Anal. Appl. 325 (2007) 542-553.
[2] A.A.M. Arafa, S.Z. Rida, M. Khalil, Fractional modeling dynamics of HIV and CD4+ T-cells during primary infection, Nonlinear Biomed. Phys. 6 (2012) 1-7.
[3] A.A.M. Arafa, S.Z. Rida, M. Khalil, The effect of anti-viral drug treatment of human immunodeficiency, Appl. Math. Modell. 37 (2013) 2189-2196.
[4] K.S. Cole, Electric conductance of biological systems, in: Proc. Cold Spring Harbor Symp. Quant. Biol, Cold Spring Harbor, New York, 1993, pp. 107-116.
[5] R.V. Culshaw, S. Ruan, A delay-differential equation model of HIV infection of CD4+ T -cells, Math. Biosci. 165 (2000) 27-39.
[6] Y. Ding, Haiping Ye, A fractional-order differential equation model of HIV infection of CD C T-cells, Math. Comput. Modell. 50 (2009) 386-392.
[7] V.D. Djordjevic, J. Jaric, B. Fabry, J.J. Fredberg, D. Stamenovic, Fractional derivatives embody essential features of cell rheological behavior, Ann. Biomed. Eng. 31 (2003) 692699.
[8] A.M.A. El-Sayed, S.Z. Rida, A.A.M. Arafa, On the solutions of time-fractional bacterial chemotaxis in a diffusion gradient chamber, Int. J. Nonlinear Sci. 7 (2009) 485-492.
[9] A.M.A. El-Sayed, A.E.M. El-Mesiry, H.A.A. El-Saka, Numerical solution for multi-term fractional (arbitrary) orders differential equations, Comput. Appl. Math. 23 (2004) 33-54.
[10] I. Hashim, O. Abdulaziz, S. Momani, Homotopy analysis method for fractional IVPs, Commun. Nonlinear Sci. Numer. Simul. 14 (2009) 674-684.
[11] I.S. Jesus, J.A.T. Machado, J.B. Cunha, Fractional electrical impedances in botanical elements, J. Vib. Control 14 (2008) 1389-1402.
[12] I.S. Jesus, J.A.T. Machado, J.B. Cunha, Fractional order electrical impedance of fruits and vegetables, in: Proceedings of the 25th IASTED International Conference Modeling, Identification, and Control, February 6-8, 2006, Lanzarote, Canary Islands, Spain, 2006.
[13] W. Lin, Global existence theory and chaos control of fractional differential equations, J. Math. Anal. Appl. 332 (2007) 709-726.
[14] M. Merdan, Tahir Khan Yev, Homotopy perturbation method for solving viral dynamical model, C.U. Fen-EdebiyatFakultesi, Fen BilimleriDergisi 31 (2010) 65-77.
[15] P.W. Nelson, A.S. Perelson, Mathematical analysis of delay differential equation models of HIV-1 infection, Math. Biosci. 179 (2002) 73.
[16] Z. Odibat, Shaher Moamni, An algorithm for the numerical solution of differential equations of fractional order, J. Appl. Math. Inform. 26 (2008) 15-27.
[17] Z. Odibat, N. Shawagfeh, Generalized Taylor's formula, Appl. Math. Comput. 186 (2007) 286-293.
[18] A.S. Perelson, Modeling the interaction of the immune system with HIV, in: C. Castillo-Chavez (Ed.), Mathematical and Statistical Approaches to AIDS Epidemiology, Lecture Notes in Biomathematics, vol. 83, Springer, New York, 1989, p. 350.
[19] A.S. Perelson, D.E. Kirschner, R. De Boer, Dynamics of HIV infection of CD4+ T cells, Math. Biosci. 114 (1993) 81-125.
[20] L.M. Petrovic, D.T. Spasic, T.M. Atanackovic, On a mathematical model of a human root dentin, Dental Mater. 21 (2005) 125-128.
[21] S.Z. Rida, H.M. El-Sherbiny, A.A.M. Arafa, On the solution of the fractional nonlinear Schrodinger equation, Phys. Lett. A 372 (2008) 553-558.
[22] F.A. Rihan, Numerical modeling of fractional-order biological systems, Abstract Appl. Anal. 2013 (2013).
[23] L. Rong, M.A. Gilchrist, Z. Feng, A.S. Perelson, Modeling within host HIV-1 dynamics and the evolution of drug resistance: tradeoffs between viral enzyme function and drug susceptibility, J. Theoret. Biol. 247 (2007) 804-818.
[24] L. Ronga, Zhilan Fenga, Alan S. Perelsonb, Emergence of HIV-1 drug resistance during antiretroviral treatment, Bull. Math. Biol. 69 (2007) 2027-2060.
[25] P.K. Srivastava, M. Banerjee, Peeyush Chandra, Modeling the drug therapy for HIV infection, J. Biol. Syst. 17 (2009) 213-223.
[26] H.C. Tuckwell, Frederic Y.M. Wan, On the behavior of solutions in viral dynamical models, Bio Systems 73 (2004) 157-161.
[27] L. Wang, Michael Y. Li, Mathematical analysis of the global dynamics of a model for HIV infection of CD4+ T cells, Math. Biosci. 200 (2006) 44-57.