Hindawi Publishing Corporation Abstract and Applied Analysis Volume 2013, Article ID 816803,11 pages http://dx.doi.org/10.1155/2013/816803

Research Article

Numerical Modeling of Fractional-Order Biological Systems

Fathalla A. Rihan1,2

1 Department of Mathematical Sciences, College of Science, UAE University, Al Ain 15551, UAE

2 Department of Mathematics, Faculty of Science, Helwan University, Cairo 11795, Egypt

Correspondence should be addressed to Fathalla A. Rihan; frihan@uaeu.ac.ae Received 17 May 2013; Accepted 23 June 2013 Academic Editor: Ali H. Bhrawy

Copyright © 2013 Fathalla A. Rihan. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

We provide a class of fractional-order differential models of biological systems with memory, such as dynamics of tumor-immune system and dynamics of HIV infection of CD4+ T cells. Stability and nonstability conditions for disease-free equilibrium and positive equilibria are obtained in terms of a threshold parameter R0 (minimum infection parameter) for each model. We provide unconditionally stable method, using the Caputo fractional derivative of order a and implicit Euler's approximation, to find a numerical solution of the resulting systems. The numerical simulations confirm the advantages of the numerical technique and using fractional-order differential models in biological systems over the differential equations with integer order. The results may give insight to infectious disease specialists.

1. Introduction

Mathematical models, using ordinary differential equations with integer order, have been proven valuable in understanding the dynamics of biological systems. However, the behavior of most biological systems has memory or aftereffects. The modelling of these systems by fractional-order differential equations has more advantages than classical integer-order mathematical modeling, in which such effects are neglected. Accordingly, the subject of fractional calculus (i.e., calculus of integral and derivatives of arbitrary order) has gained popularity and importance, mainly due to its demonstrated applications in numerous diverse and widespread fields of science and engineering. For example, fractional calculus has been successfully applied to system biology [1-5], physics [6-9], chemistry and biochemistry [10], hydrology [11, 12], medicine [13, 14], and finance [15]. In some situations, the fractional-order differential equations (FODEs) models seem more consistent with the real phenomena than the integerorder models. This is due to the fact that fractional derivatives and integrals enable the description of the memory and hereditary properties inherent in various materials and processes. Hence there is a growing need to study and use the fractional-order differential and integral equations.

Fractional-order differential equations are naturally related to systems with memory which exists in most

biological systems. Also, they are closely related to fractals [16], which are abundant in biological systems. It has been deduced in [3] that the membranes of cells of biological organism have fractional-order electrical conductance and then are classified in groups of noninteger-order models. Fractional derivatives embody essential features of cell rheological behavior and have enjoyed greatest success in the field of rheology [17]. In this paper, we propose systems of FODEs for modeling the interactions of tumor-immune system (see Section 2) and HIV infection of CD4+ T cells with immune system (see Section 3).

However, analytical and closed solutions of these types of fractional equations cannot generally be obtained. As a consequence, approximate and numerical techniques are playing important role in identifying the solution behavior of such fractional equations and exploring their applications [18-21]. The Adomian decomposition method [22, 23], extrapolation method [24], multistep method [25], monotone iterative technique [26], and predictor-corrector approach [18, 27] are proposed to provide numerical solutions for linear and nonlinear FODEs. The monograph of [28] investigates the interconnection between fractional differential equations and classical differential equations of integer-order derivatives. In this paper, we provide a numerical approach for the resulting system using implicit Euler's scheme (see Section 4). We also study the stability and convergence of the proposed method.

We first give the definition of fractional-order integration and fractional-order differentiation [29,30]. There are several approaches to the generalization of the notion of differentiation to fractional orders, for example, Riemann-Liouville, Caputo, and Generalized Functions approaches. Let L1 = L1 [a, b] be the class of Lebesgue integrable functions on [a, b], a < b < >x>.

Definition 1. The fractional integral (or the Riemann-Liouville integral) of order ß e R+ of the function f(t), t > 0

(f: R+

is defined by

I"f{t) = W)l {t-s)Mf{s)ds> t>0- (1)

The fractional derivative of order a e (n - 1,n) of f(t) is defined by two (nonequivalent) ways:

(i) Riemann-Liouville fractional derivative: take fractional integral of order (n - a), and then take nth derivative as follows:

DaJ(t) = DnJn-a fit), D: = ^n, n=l,2,.

(ii) Caputo-fractional derivative: take nth derivative, and then take a fractional integral of order (n - a)

D"f(t) = irD:f(t), n= 1,2, —

We notice that the definition of time-fractional derivative of a function f(t) at t = tn involves an integration and calculating time-fractional derivative that requires all the past history, that is, all the values of f(t) from t = 0 to t = tn. 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. The following remark addresses some of the main properties of the fractional derivatives and integrals (see [8, 29-31]).

Remark 2. Let ¡3,y e R+ and a e (0,1). Then

(i) if ll : L1 ^ L1 and f(t) e L1, then I^^f(t) =

tyf«);

(ii) lim^n iPf(x) = inf(t) uniformly on [a, b], n= 1,2, 3,...,whereI1J(t) = fQ f(s)ds;

(iii) lim^QlPfW = f(t) weakly;

(iv) if f(t) is absolutely continuous on [a,b], then lima^iD:f(t) = df(t)/dt;

(v) thus D"f(t) = (d/dt)ll-af(t) (Riemann-Liouville sense) and Daf(t) = l\-a(d/dt)f(t) (Caputo sense).

The generalized mean value theorem and another property are defined in the following remark [32].

Remark 3.

(i) Suppose that f(t) e C[a,b] and DaJ(t) e C(a,b] for 0 < a < 1, and then we have

f(t) = f(a) + -^D:f(V(t-aT,

r (a) (4)

with a < £ < t Vt e (a,b].

(ii) If (i) holds and DaJ(t) > 0 for all t e [a,b], then f(t) is nondecreasing for each t e [a, b]. If D"f(t) < 0 for all t e [a, b], then f(t) is nonincreasing for each t e [a,b].

We next provide a class of fractional-order differential models to describe the dynamics of tumour-immune system interactions.

2. Fractional Model of Tumor-Immune System

Immune system is one of the most fascinating schemes from the point of view of biology and mathematics. The immune system is complex, intricate, and interesting. It is known to be multifunctional and multipathway, so most immune effectors do more than one job. Also each function of the immune system is typically done by more than one effector, which makes it more robust. The reason of using FODEs is that they are naturally related to systems with memory which exists in tumor-immune interactions.

Ordinary and delay differential equations have long been used in modeling cancer phenomena [33-37], but fractional-order differential equations have short history in modeling such systems with memory. The authors in [1] used a system of fractional-order differential equations in modeling cancer-immune system interaction. The model includes two immune effectors: E1(t), E2(t) (such as cytotoxic T cells and natural killer cells), interacting with the cancer cells, T(t), with a Holling function of type III. (Holling type III describes a situation in which the number of prey consumed per predator initially rises slowly as the density of prey increases but then levels off with further increase in prey density. In other words the response of predators to prey is depressed at low prey density, then levels off with further increase in prey density.) The model takes the form

DaT = aT-r1TE1 - r2TE2,

DaE1 =

DaE, =

-diEi + T2 + ik , 0 < a < 1,

-d2E2 + ■ r , 22 T2 +k1

where T = T(t), E1 = E1(t), E2 = E2(t), and a, r1, r2, d1, d2, k1, and k2 are positive constants. The interaction terms in the second and third equations of model (5) satisfy the crossreactivity property of the immune system. It has been assumed that (d1k1/(1 - d1)) < (d2k2/(1 - d2)) to avoid

the nonbiological interior solution where both immune effectors coexist. The equilibrium points of the system (5) are

E0 = (0, 0,0), E1 =

i1 k1 a

(l-d1)'r1

2-2 Q a

{l-d2)' ' Г2

The first equilibrium E0 is the nave, the second E1 is the memory, and the third E2 is endemic according to the value of the tumor size. Stability analysis shows that the nave state is unstable. However, the memory state is locally asymptotically stable if dx < d2 and dx < 1. While the endemic state is locally asymptotically stable if d2 < d1 and d2 < 1, there is bifurcation at dx = 1. The stability of the memory state depends on the value of one parameter, namely, the immune effector death rate.

Now we modify model (5) to include three populations of the activated immune-system cells, E(t); the tumor cells, T(t); and the concentration of IL-2 in the single tumor-site compartment, IL(t). We consider the classic bilinear model that includes Holling function of type I and external effector cells ^ and input of IL-2, s2. (holling Type I is a linear relationship, where the predator is able to keep up with increasing density of prey by eating them in direct proportion to their abundance in the environment. If they eat 10% of the prey at low density, they continue to eat 10% of them at high densities.) The interactions of the three populations are then governed by the fractional-order differential model:

DXl E = s1 + p1ET - p2E + p3EIL,

D"2T = p4T (1 - p5T) - p6ET, 0 < at <1, i = 1,2,3,

D"3 IL = S2 + p7ET-ps IL,

with initial conditions: E(0) = E0, T(0) = T0, and IL(0) = ILo. The parameter p1 is the antigenicity rate of the tumor (immune response to the appearance of the tumor), and ^ is the external source of the effector cells, with rate of death p2. The parameter p4 incorporates both multiplication and death of tumor cells. The maximal carrying capacity of the biological environment for tumor cell is P5-1; p3 is considered as the cooperation rate of effector cells with interleukin-2 parameter; p6 is the rate of tumor cells; p7 is the competition rate between the effector cells and the tumor cells. External input of IL-2 into the system is s2, and the rate loss parameter of IL-2 cells is p8.

In the absence of immunotherapy with IL-2, we have

D"1 E = s1 +p1 ET-p2E,

0 <ai <1, i = 1,2, (8)

Da2 T = p4T(l-p5T)-p6ET.

To ease the analysis and stability of the steady states with meaningful parameters and minimize sensitivity (or robustness) ofthe model, we nondimensionalize the bilinear system (8) by taking the rescaling

Eo' P2

pJo tsEo '

ь = Ps To,

Therefore, after the previous substitution into (8) and replacing t by t, the model becomes

D"1 x = a + шху - вх, D"2у = ay(\- by) - xy.

2.1. Equilibria and Local Stability of Model (10). The steady states of the reduced model (10) are again the intersection of the null clines D"1 x = 0, D"2у = 0. If у = 0, the tumorfree equilibrium is at E0 = (x,~y) = (а/в, 0). This steady state is always exist, since а/в > 0. From the analysis, it is easy to prove that the tumor-free equilibrium E0 = (а/в, 0) of the model (10) is asymptotically stable if threshold parameter (the minimum tumor-clearance parameter) R0 = ав/а < 1 and unstable if R0 >1.

However, if у = 0, the steady states are obtained by solving

waby - a(w + вЪ) у + ав - a = 0.

(7) In this case, we have two endemic equilibria, E1 and E2:

E1 = (x1,y1), where x1 =

En = (х2,Уп), where ^ =

-а (Ъв -w)- VÄ 2w

а(Ьв + w) + VÄ 2abw

-а (Ьв -w) + VÄ 2w

а(Ьв + w) - VÄ 2abw

with A = a2(bd - w)2 + Aawab > 0. The Jacobian matrix of the system (10) at the endemic equilibrium E1 is

j(E1) = (iüy1--e

V -y1 a-2aby1 -x1J

T = tst

Proposition 4. Assume that the endemic equilibrium E1 exists and has nonnegative coordinates. If R0 = ad/a < 1, then tr(J(E1)) > 0 and E1 is unstable.

Proof. Since

tr (J(Ei)) =

o2 -o (ab + bd) - ab2d

o - ab 2abo

ya2(bd + o)2 - 4abo (ad - a),

then inequality tr(/(E1)) > 0 is true if a [o2 - co (ab + bd) - ab2d]

> (ab - o) ^a2(bd + o)2 - 4abo (ad - a).

Therefore, when ad < <r,wehave co2- ob(a+d)-adb2 > 0,and hence both sides of the inequality are positive. Therefore if the equilibrium point E1 exists and has nonnegative coordinates, then tr(J(E1)) > 0 and the point (E1) is unstable, whenever

R0 = ad/a < 1.

Similarly, we arrive at the following proposition.

Proposition 5. If the point E2 exists and has nonnegative coordinates, then it is asymptotically stable.

We next examine a fractional-order differential model for HIV infection of CD4+ T cells.

3. Fractional Model of HIV Infection of CD4+ T Cells

As mentioned before, many mathematical models have been developed to describe the immunological response to infection with human immunodeficiency virus (HIV), using ordinary differential models (see, e.g., [38, 39]) and delay differential models [40, 41]. In this section, we extend the analysis and replace the original system of Perelson et al. [39] by an equivalent fractional-order differential model of HIV infection of CD4+ T cells that consists of three equations:

DaiH = s - pHH + rH ( 1 -

) - k1VH,

D"2I = k[VH - fal, Da,V = MphI - k1VH -

0.5 <(X: <1, i = 1,2,3.

Here H = H(t) represents the concentration of healthy CD4+ T cells at time t, I = I(t) represents the concentration of infected CD4+ T cells, and V = V(t) is the concentration of free HIV at time t. In this system, the logistic growth of

the healthy CD4+ T-cells is (1 - ((H + I)/Hmax)), and the proliferation of infected cells is neglected. The parameter s is the source of CD4+ T cells from precursors, is the natural death rate of CD4+ T cells (^HHmax > s, cf. [39, page 85]), r is their growth rate (thus, r > ^H in general), and Hmax is their carrying capacity. The parameter k1 represents the rate of infection of T cells with free virus. k'1 is the rate at which infected cells become actively infected. ^ is a blanket death term for infected cells to reflect the assumption that we do not initially know whether the cells die naturally or by bursting. In addition, is the lytic death rate for infected cells. Since M viral particles are released by each lysing cell, this term is multiplied by the parameter M to represent the source for free virus (assuming a one-time initial infection). Finally, is the loss rate of virus. The initial conditions for infection by free virus are H(0) = HQ, 1(0) = IQ, and V(0) = VQ.

Theorem 6. According to Remark 3 and the fact that s > 0, then the system (16) has a unique solution (H, I, V) which remains in R+ and bounded by Hmax; [11].

3.1. Equilibria and Local Stability of Model (16). To evaluate the equilibrium points of system (16), we put DaiH(t) = D"2I(t) = D"3V(t) = 0. Then the infection-free equilibrium point (uninfected steady state) is Eq = (HQ, 0,0), and endemic equilibrium point (infected steady state) is E+ = (H*,I*,V*),where

= k\H*V* Pi

V, = Pi [(s +(r- Ph) H')Hmax - rH*

H [k[rW -kMHmJ

The Jacobian matrix J(EQ) for system (16) evaluated at the uninfected steady state EQ is then given by

-pH + r - 2rH0 -rH0

- k1 H0

0 M^h - (kH + yy)

Let us introduce the following definition and assumption to ease the analysis.

Definition 7. The threshold parameter RQ (the minimum infection-free parameter) is the parameter that has the property that if RQ < 1, then the endemic infected state does not exist, while if RQ > 1, the endemic infected state persists, where

k[phMHo

Ro Pi + k1Hoï

We also assume that

2rHmax

The uninfected steady state is asymptotically stable if all of the eigenvalues X of the Jacobian matrix /(EQ ),given by (18), have negative real parts. The characteristic equation det(/(E0) -I) = 0 becomes

(x + ^H-r + ) (X2 + BX + C) = 0,

\ ^max '

where B = fa1 + k1H0 + fav and C = fa1(k1H0 + fav)-k[fahMH0. Hence, the three roots of the characteristic equation (21) are

= -Vh + r- = < 0

X23 = 1[-B±^B2-4C].

Proposition 8. IfR** = (k[fahMH0)/fa(fav + k1H0) < 1, then C > 0 and the three roots of the characteristic equation (21) will have negative real parts.

Corollary 9. In case of uninfected steady state EQ, one has three cases.

(i) If R** < 0, the uninfected state is asymptotically stable and the infected steady state E+ does not exist (unphysical).

(ii) If RQ = 1, then C = 0, andfrom (21) implies that one eigenvalue X = 0 and the remaining two eigenvalues have negative real parts. The uninfected and infected steady states collide, and there is a transcritical bifurcation.

(iii) If RQ > 1, then C < 0, and thus at least one eigenvalue will be positive real root. Thus, the uninfected state E0 is unstable, and the endemically infected state E+ emerges.

To study the local stability of the positive infected steady states EQ for RQ > 1, we consider the linearized system of (16) at E*. The Jacobian matrix at E+ becomes

j(E+) =

TT ^max

k1V* Mfa

Mpb -(kiH*+^v)

^H - r + k1V* +

r(2H* +T)

P (X) = X3 + a1X2 + a2X + a3 = 0,

Then the characteristic equation of the linearized system is

a1 = fa + fav + k1H* + V, a2 = L* (fa + fav + k1H*) + fa (fav + k1H*)

-k\H*V* - k[ H* (Mfa -

a? = k[ H*

kiMfaV* -TàL-L*Mpb

+ L*fa1 (fav + k1H*) - fa1k1H*V*.

The infected steady state EQ is asymptotically stable if all of the eigenvalues have negative real parts. This occurs if and only if Routh-Hurwitz conditions are satisfied; that is, a1 > 0, a3 > 0, and a1a2 > a3.

4. Implicit Euler's Scheme for FODEs

Since most of the fractional-order differential equations do not have exact analytic solutions, approximation and numerical techniques must be used. Several numerical methods have been proposed to solve the fractional-order differential equations [18,42,43]. In addition, most of resulting biological systems are stiff. (One definition of the stiffness is that the global accuracy of the numerical solution is determined by stability rather than local error, and implicit methods are more appropriate for it.) The stiffness often appears due to the differences in speed between the fastest and slowest components of the solutions and stability constraints. In addition, the state variables of these types of models are very sensitive to small perturbations (or changes) in the parameters occur in the model. Therefore, efficient use of a reliable numerical method that based in general on implicit formulae for dealing with stiff problems is necessary.

Consider biological models in the form of a system of FODEs of the form

DaX (t) = F(t,X(t), P), i e [0, J], 0 < a < 1, X(0) = Xq.

Here X(t) = [x1(t),x2(t),..., xn(t)]T, P is the set of parameters appear in the model, and F(t, X(t)) satisfies the Lipschitz condition

\\F (t, X(t), P)-F (t, Y (t), P)|| < ^yX(i) - Y (t)y,

where Y(t) is the solution of the perturbed system.

Theorem 10. Problem (27) has a unique solution provided that the Lipschitz condition (28) is satisfied and KTa/T(a+ 1) < 1.

Proof. Using the definitions of Section 1, we can apply a fractional integral operator to the differential equation (27)

and incorporate the initial conditions, thus converting the equation into the equivalent equation

X(t) = X (0) + f (t - s)a-1F (s, X (s), P) ds (29) r (a) Jq

which also is a Volterra equation of the second kind. Define the operator L, such that

LX (t) = X(0) + -^ f (t - s)a~LF (s, X (s), P) ds. r (a) Jo

\a-1 i

Then, we have \\LX(t)- LY(t)\\

— I (t-sf-1 \\F(s,X(s), P) - F (s, Y (s), (a) o

\ (t-s)a-1 sup \X(s)-Y(s)\ds Jo seio.ri

r (a) Jo K "

P)\\ds

r(a) J0 se[0,Tl

<-^\\X-Y\\ I sa-1ds r (a) Jo

<-^\\X-Y\\Ta. r (a)

Thus for KTa/T(a + 1) < 1, we have

\\LX(t)- LY(t)\\<\\X-Y\\. This implies that our problem has a unique solution.

However, converted Volterra integral equation (29) is with a weakly singular kernel, such that a regularization is not necessary any more. It seems that there exists only a very small number of software packages for nonlinear Volterra equations. In our case the kernel may not be continuous, and therefore the classical numerical algorithms for the integral part of (29) are unable to handle the solution of (27). Therefore, we implement the implicit Euler's scheme to approximate the fractional-order derivative.

Given model (27) and mesh points T = |iQ, t1,..., tN], such that tQ = 0 and tN = T. Then a discrete approximation to the fractional derivative can be obtained by a simple quadrature formula, using the Caputo fractional derivative (3) of order a, 0 < a < 1, and using implicit Euler's approximation as follows (see [44]):

D>, (tn)

1 I dxj (s) r (1 - a) Jo ds

(tn - s) adü

1_I (jh

-a)p[)(j-i

xj -xj-1

+ 0(h)

(nh - s) ads

(1 -a)r(1 - a)

X -xj-1

+ 0(h)

[{n-j+1)1 * -(n-j)1 a]

(1 -a)r(1 - a) ha

-I[xj-xf1] [(n- j +1)1--(n- j)1-] j=1

(1 -a)r(1 -a)

-I [xj - xj-1 ][(n-j+ 1)1-a - (n - j)1-"] 0 (h2-*) . j=1

Setting G (a, h) =

(1 -a)r(1 - a) ha

■ 1-a / ■ , \1-a = j -(j-1) ,

(where o* = 1), (34)

then the first-order approximation method for the computation of Caputo's fractional derivative is then given by the expression

Dy, (tn) = G (a, h) Y"* (*7i+1 - + O(h). (35)

From the analysis and numerical approximation, we also arrive at the following proposition.

Proposition 11. The presence of a fractional differential order in a differential equation can lead to a notable increase in the complexity of the observed behaviour, and the solution is continuously depends on all the previous states.

4.1. Stability and Convergence. We here prove that the fractional-order implicit difference approximation (35) is unconditionally stable. It follows then that the numerical solution converges to the exact solution as h — 0.

In order to study the stability of the numerical method, let us consider a test problem of linear scaler fractional differential equation

DQu(t) = pQu(t) + p1, U(0) = Uq, (36)

such that 0 < a < 1, and pQ < 0, p1 > 0 are constants.

È^ 0.5

ft \ hit) a = 0.75

I T(t)

20 40 60

Bq" 1.5

a = 0.75

Ï\ El (t)

I T(t)

)(t 1.2

Kf 0.8

a = 0.95

E2 (t)

1.6 1.4 1.2 1

0.8 0.6 0.4 0.2 0

a = 0.75

Figure 1: Numerical simulations of the FODEs model (5) when a = 0.75 and a = 0.95 with (a = r1 = r2 = 1; d1 = 0.3; d2 = 0.7; k1 = 0.3; k2 = 0.7) andwhen dx = 0.7; d2 = 0.3. The system converges to the stable steady states Ep E2. The fractional derivative damps the oscillation behavior.

Theorem 12. The fully implicit numerical approximation (35), Since (1 - (p0/Gah)) > 1 for all Gah, then to test problem (36) for all t > 0, is consistent and unconditionally stable.

Proof. We assume that the approximate solution of (36) is of the form u(tn) ~Un = ^n, and then (36) can be reduced to

1--^ R

= t*-i + I«'f (^J -tn-j+i) + Gr, n^2

„ t.,-1 + IU "['f (tn-j - L-j+1 ) + Pi IG «h

C,n = -r-, »-> n>2.

(1 - (PQIG«,))

Thus, for n = 2, the previous inequality implies

C2 + ^ «0 -Q-

C < tn-i + (!>«-} - tn-j+1 ), n>2. (40)

Using the relation (39) and the positivity of the coefficients w2, we get

Repeating the process, we have from (40)

Cn < tn-i + (tn-j - tn-j+1 )<tn-i> (43)

4 3.5 3 2.5 2 1.5 1

a = 1 a = 0.9 a = 0.7

b A A „ .

1 b ' \ [ \J V v w

a = 1 a = 0.9 a = 0.7

a = 1 a = 0.9 a = 0.7

1.5 2 2.5 3 Effector cells, E(t)

Figure 2: Numerical simulations of the FODEs model (10) when s = 0.1181, " = 0.1184, 0 = 0.1747, r = 0.636, and b = 0.002. a1 = a2 = 1,0.9,0.7. The endemic state E2 is locally asymptotically stable when RQ = a0/a > 1. The fractional derivative damps the oscillation behavior.

since each term in the summation is negative. Thus Cn < Cn-1 < (n-2 < ■■■ < CQ. With the assumption that (n = \Un\ < cQ = \UQ\, which entails \\Un\\ < \\UQ\, we have stability. □

Of course this numerical technique can be used both for linear and for nonlinear problems, and it may be extended to multiterm FODEs. For more details about stability and convergence of the fractional Euler method, we refer to [19, 45].

4.2. Numerical Simulations. We employed the implicit Euler's scheme (35) to solve the resulting biological systems of FODEs (5), (10), and (16). Interesting numerical simulations of the fractional tumor-immune models (5) and (10), with step size h = 0.05 and 0.5 < a < 1 and parameters values given in the captions, are displayed in Figures 1, 2, and 3. The

equilibrium points (for infection-free and endemic cases) are the same in both integer-order (when a = 1) and fractional-order (a < 1) models. We notice that, in the endemic steady states, the fractional-order derivative damps the oscillation behavior. From the graphs, we can see that FODEs have rich dynamics and are better descriptors of biological systems than traditional integer-order models.

The numerical simulations for the HIV FODE model (16) (with parameter values given in the captions) are displayed in Figure 4. We note that the solution of the model, with various values of a, continuously depends on the time-fractional derivative but arrives to the equilibrium points. The displayed solutions, in the figure, confirm that the fractional order of the derivative plays the role of time delay in the system. We should also note that although the equilibrium points are the same for both integer-order and fractional-order models,

0 20 40 60 80 100 120 140 160 180 200 Time

— a = 1

■ a = 0.9

■ a = 0.8

0 20 40 60 80 100 120 140 160 180 200 Time

- a = 1

■ a = 0.9

■ a = 0.8

0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 Effector cells, E(t)

- a = 1

■ a = 0.9

■ a = 0.8

Figure 3: Numerical simulations of the FODEs model (10) when s = 0.1181, w = 0.1184, 0 = 0.3747, r = 1.636, and b = 0.002. al = a2 1,0.9,0.7. The infection-free steady state E0 is locally asymptotically stable when R0 = a0/a < 1.

the solution of the fractional order model tends to the fixed point over a longer period of time.

5. Conclusions

In fact, fractional-order differential equations are generalizations of integer-order differential equations. Using fractional-order differential equations can help us to reduce the errors arising from the neglected parameters in modeling biological systems with memory and systems distributed parameters. In this paper, we presented a class of fractional-order differential models of biological systems with memory to model the interaction of immune system with tumor cells and with HIV infection of CD4+ T cells. The models possess nonnegative solutions, as desired in any population dynamics. We obtained the threshold parameter R0 that represents the minimum tumor-clearance parameter or minimum infection free for each model.

We provided unconditionally stable numerical technique, using the Caputo fractional derivative of order a and implicit Euler's approximation, for the resulting system. The numerical technique is suitable for stiff problems. The solution of the system at any time t* is continuously depends on all the previous states at t < t*. We have seen that the presence of the fractional differential order leads to a notable increase in the complexity of the observed behaviour and play the role of time lag (or delay term) in ordinary differential model. We have obtained stability conditions for disease-free equilibrium and nonstability conditions for positive equilibria. We should also mention that one of the basic reasons of using fractional-order differential equations is that fractional-order differential equations are, at least, as stable as their integer-order counterpart. In addition, the presence of a fractional differential order in a differential equation can lead to a notable increase in the complexity of the observed behaviour, and the solution continuously depends on all

1200 .§ 1000

lu £ 200

a = 0.9

[a = 0.7 ;

0 20 40 60 80

100 120 140 160 180 200 Time

x10 10

o £ c

\a = 0.9

/ ; J\ ; ; ;

0 20 40 60 80

100 120 Time

140 160 180 200

T3 500

a 300 Q

V\ a = 0.9 ;

IV = 0.7 - - V - -

j\ ; J\

0 20 40 60 80

100 120 Time

140 160 180 200

400 600 H(t)

Figure 4: Numerical simulations of the FODEs model (16) when s = 20 day 1 mm 3, = 0.02 day 1, = 0.26 day 1, = 0.24 day 1, = 2.4 day-1, k1 = 2.4 x 10-5 mm3 day-1, k'l = 2x 10-5 mm3 day-1, r = 0.3 day-1, Hmax = 1500 mm-3, and M = 1400. The trajectories of the system approach to the steady state E**.

the previous states. The analysis can be extended to other

models related to the immune response systems.

Acknowledgments

The work is funded by the UAE University, NRF Project no.

NRF-7-20886. The author is grateful to Professor Changpin

Li's valuable suggestions and referees' constructive comments.

References

[1] E. Ahmed, A. Hashish, and F. A. Rihan, "On fractional order cancer model," Journal of Fractional Calculus and Applied Analysis, vol. 3, no. 2, pp. 1-6, 2012.

[2] A. A. M. Arafa, S. Z. Rida, and M. Khalil, "Fractional modeling dynamics of HIV and CD4+ T-cells during primary infection," Nonlinear Biomedical Physics, vol. 6, no. 1, article 1, 2012.

[3] K. S. Cole, "Electric conductance of biological systems," Cold Spring Harbor Symposia on Quantitative Biology, pp. 107-116, 1993.

[4] A. M. A. El-Sayed, A. E. M. El-Mesiry, and H. A. A. El-Saka, "On the fractional-order logistic equation," Applied Mathematics Letters, vol. 20, no. 7, pp. 817-823, 2007.

[5] H. Xu, "Analytical approximations for a population growth model with fractional order," Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 5, pp. 1978-1983, 2009.

[6] L. Debnath, "Recent applications of fractional calculus to science and engineering," International Journal of Mathematics and Mathematical Sciences, no. 54, pp. 3413-3442, 2003.

[7] A. M. A. El-Sayed, "Nonlinear functional-differential equations of arbitrary orders," Nonlinear Analysis. Theory, Methods & Applications, vol. 33, no. 2, pp. 181-186, 1998.

[8] R. Hilfer, Ed., Applications of Fractional Calculus in Physics, World Scientific, River Edge, NJ, USA, 2000.

[9] G. M. Zaslavsky, "Chaos, fractional kinetics, and anomalous transport," Physics Reports, vol. 371, no. 6, pp. 461-580, 2002.

[10] S. B. Yuste, L. Acedo, and K. Lindenberg, "Reaction front in an A + B ^ C reaction-subdiffusion process," Physical Review E, vol. 69, no. 3, Article ID 036126, pp. 1-36126, 2004.

[11] W. Lin, "Global existence theory and chaos control of fractional differential equations," Journal of Mathematical Analysis and Applications, vol. 332, no. 1, pp. 709-726, 2007.

[12] H. Sheng, Y. Q. Chen, and T. S. Qiu, Fractional Processes and Fractional-Order Signal Processing, Springer, New York, NY, USA, 2012.

[13] K. Assaleh and W. M. Ahmad, "Modeling of speech signals using fractional calculus," in Proceedings of the 9th International Symposium on Signal ProcessinganditsApplications (ISSPA '07), Sharjah, United Arab Emirates, February 2007.

[14] Y. Ferdi, "Some applications of fractional order calculus to design digital filters for biomedical signal processing," Journal of Mechanics in Medicine and Biology, vol. 12, no. 2, Article ID 12400088, 13 pages, 2012.

[15] W.-C. Chen, "Nonlinear dynamics and chaos in a fractional-order financial system," Chaos, Solitons and Fractals, vol. 36, no. 5, pp. 1305-1314, 2008.

[16] A. Rocco and B. J. West, "Fractional calculus and the evolution of fractal phenomena," Physica A, vol. 265, no. 3, pp. 535-546, 1999.

[17] V. D. Djordjevic, J. Jaric, B. Fabry, J. J. Fredberg, and D. Stamenovic, "Fractional derivatives embody essential features of cell rheological behavior," Annals of Biomedical Engineering, vol. 31, no. 6, pp. 692-699, 2003.

[18] K. Diethelm, N. J. Ford, and A. D. Freed, "A predictor-corrector approach for the numerical solution of fractional differential equations," Nonlinear Dynamics, vol. 29, no. 1-4, pp. 3-22,2002.

[19] C. Li and F. Zeng, "The finite difference methods for fractional ordinary differential equations," Numerical Functional Analysis and Optimization, vol. 34, no. 2, pp. 149-179, 2013.

[20] J. C. Trigeassou, N. Maamri, J. Sabatier, and A. Oustaloup, "State variables and transients of fractional order differential systems," Computers & Mathematics with Applications, vol. 64, no. 10, pp. 3117-3140, 2012.

[21] X. Zhang and Y. Han, "Existence and uniqueness of positive solutions for higher order nonlocal fractional differential equations," Applied Mathematics Letters, vol. 25, no. 3, pp. 555-560, 2012.

[22] A. M. A. El-Sayed and M. Gaber, "The Adomian decomposition method for solving partial differential equations of fractal order in finite domains," Physics Letters A, vol. 359, no. 3, pp. 175-182, 2006.

[23] H. Jafari and V. Daftardar-Gejji, "Solving a system of nonlinear fractional differential equations using Adomian decomposition," Journal of Computational and Applied Mathematics, vol. 196, no. 2, pp. 644-651, 2006.

[24] K. Diethelm and G. Walz, "Numerical solution of fractional order differential equations by extrapolation," Numerical Algorithms, vol. 16, no. 3-4, pp. 231-253,1997.

[25] L. Galeone and R. Garrappa, "On multistep methods for differential equations of fractional order," Mediterranean Journal of Mathematics, vol. 3, no. 3-4, pp. 565-580, 2006.

[26] G. Wang, R. P. Agarwal, and A. Cabada, "Existence results and the monotone iterative technique for systems of nonlinear fractional differential equations," Applied Mathematics Letters, vol. 25, no. 6, pp. 1019-1024, 2012.

[27] N. J. Ford and A. C. Simpson, "The numerical solution of fractional differential equations: speed versus accuracy," Numerical Algorithms, vol. 26, no. 4, pp. 333-346, 2001.

[28] V. Lakshmikantham, S. Leela, and J. Vasundhara, Theory of Fractional Dynamic Systems, Cambridge Academic Publishers, Cambridge. UK.

[29] I. Podlubny, Fractional Differential Equations, vol. 198 of Mathematics in Science and Engineering, Academic Press, San Diego, Calif, USA, 1999.

[30] S. G. Samko, A. A. Kilbas, and O. I. Marichev, Fractional integrals and derivatives, Gordon and Breach Science, Yverdon, Switzerland, 1993.

[31] C. Li and Y. Ma, "Fractional dynamical system and its linearization theorem," Nonlinear Dynamics, vol. 71, no. 4, pp. 621-633, 2013.

[32] Z. M. Odibat and N. T. Shawagfeh, "Generalized Taylor's formula," Applied Mathematics and Computation, vol. 186, no. 1, pp. 286-293, 2007.

[33] N. Bellomo, A. Bellouquid, J. Nieto, and J. Soler, "Multiscale biological tissue models and flux-limited chemotaxis for multicellular growing systems," Mathematical Models & Methods in Applied Sciences, vol. 20, no. 7, pp. 1179-1207, 2010.

[34] A. Gokdogan, A. Yildirim, and M. Merdan, "Solving a fractional order model of HIV infection of CD+ T cells," Mathematical and Computer Modelling, vol. 54, no. 9-10, pp. 2132-2138, 2011.

[35] D. Kirschner and J. C. Panetta, "Modeling immunotherapy of the tumor—immune interaction," Journal of Mathematical Biology, vol. 37, no. 3, pp. 235-252,1998.

[36] F. A. Rihan, M. Safan, M. A. Abdeen, and D. Abdel Rahman, "Qualitative and computational analysis of a mathematical model for tumor-immune interactions," Journal of Applied Mathematics, vol. 2012, Article ID 475720,19 pages, 2012.

[37] R. Yafia, "Hopf bifurcation in differential equations with delay for tumor-immune system competition model," SIAM Journal on Applied Mathematics, vol. 67, no. 6, pp. 1693-1703, 2007.

[38] K. A. Pawelek, S. Liu, F. Pahlevani, and L. Rong, "A model of HIV-1 infection with two time delays: mathematical analysis and comparison with patient data," Mathematical Biosciences, vol. 235, no. 1, pp. 98-109, 2012.

[39] A. S. Perelson, D. E. Kirschner, and R. De Boer, "Dynamics of HIV infection of CD4+ T cells," Mathematical Biosciences, vol. 114, no. 1, pp. 81-125,1993.

[40] G. A. Bocharov and F. A. Rihan, "Numerical modelling in biosciences using delay differential equations," Journal of Computational and Applied Mathematics, vol. 125, no. 1-2, pp. 183199, 2000.

[41] R. V. Culshaw and S. Ruan, "A delay-differential equation model of HIV infection of CD4+ T-cells," Mathematical Biosciences, vol. 165, no. 1, pp. 27-39, 2000.

[42] R. Anguelov and J. M.-S. Lubuma, "Nonstandard finite difference method by nonlocal approximation," Mathematics and Computers in Simulation, vol. 61, no. 3-6, pp. 465-475, 2003.

[43] C. Li and F. Zeng, "Finite difference methods for fractional differential equations," International Journal of Bifurcation and Chaos in Applied Sciences and Engineering, vol. 22, no. 4, 28 pages, 2012.

[44] F. A. Rihan, "Computational methods for delay parabolic and time-fractional partial differential equations," Numerical Methods for Partial Differential Equations, vol. 26, no. 6, pp. 1556-1571, 2010.

[45] Z. Odibat and S. Momani, "Numerical methods for nonlinear partial differential equations of fractional order," Applied Mathematical Modelling, vol. 32, no. 1, pp. 28-39, 2008.

Copyright of Abstract & Applied Analysis is the property of Hindawi Publishing Corporation and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.