# Numerical Modeling of Fractional-Order Biological Systems

Authors: Fathalla A. Rihan
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

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 ℛ0 (minimum infection parameter) for each model. We provide unconditionally stable method, using the Caputo fractional derivative of order α 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.


---

## Body

## 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 integer-order 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<∞.Definition 1.
The fractional integral (or the Riemann-Liouville integral) of orderβ∈ℝ+ of the function f(t), t>0 (f:ℝ+→ℝ) is defined by
(1)Iaβf(t)=1Γ(β)∫at‍(t-s)β-1f(s)ds,t>0.
The fractional derivative of order α∈(n-1,n) of f(t) is defined by two (nonequivalent) ways: (i)
Riemann-Liouville fractional derivative: take fractional integral of order(n-α), and then take nth derivative as follows:
(2)D*αf(t)=D*nIan-αf(t),D*n=dndtn,n=1,2,….(ii)
Caputo-fractional derivative: takenth derivative, and then take a fractional integral of order (n-α)(3)Dαf(t)=Ian-αD*nf(t),n=1,2,….We notice that the definition of time-fractional derivative of a functionf(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β,γ∈ℝ+ and α∈(0,1). Then (i)
ifIaβ:L1→L1 and f(t)∈L1, then IaβIaγf(t)=Iaβ+γf(t);(ii)
limβ→nIaβf(x)=Ianf(t) uniformly on [a,b], n=1,2,3,…, where Ia1f(t)=∫0t‍f(s)ds;(iii)
limβ→0Iaβf(t)=f(t) weakly;(iv)
iff(t) is absolutely continuous on [a,b], then limα→1D*αf(t)=df(t)/dt;(v)
thusD*αf(t)=(d/dt)I*1-αf(t) (Riemann-Liouville sense) and Dαf(t)=I*1-α(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 thatf(t)∈C[a,b] and D*αf(t)∈C(a,b] for 0<α≤1, and then we have
(4)f(t)=f(a)+1Γ(α)D*αf(ξ)(t-a)α,witha<ξ<t∀t∈(a,b].
(ii)
If (i) holds andD*αf(t)≥0 for all t∈[a,b], then f(t) is nondecreasing for each t∈[a,b]. If D*αf(t)≤0 for all t∈[a,b], then f(t) is nonincreasing for each t∈[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
(5)DαT=aT-r1TE1-r2TE2,DαE1=-d1E1+T2E1T2+k1,0<α≤1,DαE2=-d2E2+T2E2T2+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
(6)ℰ0=(0,0,0),ℰ1=(d1k1(1-d1),ar1,0),ℰ2=(d2k2(1-d2),0,ar2).
The first equilibrium ℰ0 is the nave, the second ℰ1 is the memory, and the third ℰ2 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 d1<d2 and d1<1. While the endemic state is locally asymptotically stable if d2<d1 and d2<1, there is bifurcation at d1=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 s1 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:
(7)Dα1E=s1+p1ET-p2E+p3EIL,Dα2T=p4T(1-p5T)-p6ET,0<αi≤1,i=1,2,3,Dα3IL=s2+p7ET-p8IL,
with initial conditions: E(0)=E0, T(0)=T0, and IL(0)=IL0. The parameter p1 is the antigenicity rate of the tumor (immune response to the appearance of the tumor), and s1 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(8)Dα1E=s1+p1ET-p2E,Dα2T=p4T(1-p5T)-p6ETa0<αi≤1,i=1,2,Dα2T=p4T(1-p5T)-p6ET.
To ease the analysis and stability of the steady states with meaningful parameters and minimize sensitivity (or robustness) of the model, we nondimensionalize the bilinear system (8) by taking the rescaling
(9)x=EE0,y=TT0,ω=p1T0tsE0,θ=p2ts,σ=s1tsE0,a=p4ts,b=p5T0,1=p6E0ts,τ=tst.
Therefore, after the previous substitution into (8) and replacing τ by t, the model becomes
(10)Dα1x=σ+ωxy-θx,Dα2y=ay(1-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α1x=0, Dα2y=0. If y=0, the tumor-free equilibrium is at ℰ¯0=(x¯,y¯)=(σ/θ,0). This steady state is always exist, since σ/θ>0. From the analysis, it is easy to prove that the tumor-free equilibrium ℰ¯0=(σ/θ,0) of the model (10) is asymptotically stable if threshold parameter (the minimum tumor-clearance parameter) ℛ¯0=aθ/σ<1 and unstable if ℛ¯0>1.However, ify≠0, the steady states are obtained by solving
(11)ωaby2-a(ω+θb)y+aθ-σ=0.
In this case, we have two endemic equilibria, ℰ¯1 and ℰ¯2:
(12)ℰ¯1=(x¯1,y¯1),wherex¯1=-a(bθ-ω)-Δ2ω,ℰ¯1=(x¯1,y¯1),wherex¯1y¯1=a(bθ+ω)+Δ2abω,ℰ¯2=(x¯2,y¯2),wherex¯2=-a(bθ-ω)+Δ2ω,ℰ¯2=(x¯2,y¯2),wherex¯2y¯2=a(bθ+ω)-Δ2abω
with Δ=a2(bθ-ω)2+4σωab>0. The Jacobian matrix of the system (10) at the endemic equilibrium ℰ¯1 is
(13)J(ℰ¯1)=(ωy¯1-θωx¯1-y¯1a-2aby¯1-x¯1).Proposition 4.
Assume that the endemic equilibriumℰ¯1 exists and has nonnegative coordinates. If ℛ¯0=aθ/σ<1, then tr(J(ℰ¯1))>0 and ℰ¯1 is unstable.Proof.
Since(14)tr(J(ℰ¯1))=ω2-ω(ab+bθ)-ab2θ2bωtr(J(ℰ¯1))=+ω-ab2abωa2(bθ+ω)2-4abω(aθ-σ),
then inequality tr(J(ℰ¯1))>0 is true if
(15)a[ω2-ω(ab+bθ)-ab2θ]>(ab-ω)a2(bθ+ω)2-4abω(aθ-σ).
Therefore, when aθ<σ, we have ω2-ωb(a+θ)-aθb2>0, and hence both sides of the inequality are positive. Therefore if the equilibrium point ℰ¯1 exists and has nonnegative coordinates, then tr(J(ℰ¯1))>0 and the point (ℰ¯1) is unstable, whenever ℛ¯0=aθ/σ<1.Similarly, we arrive at the following proposition.Proposition 5.
If the pointℰ¯2 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.


## 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α1x=0, Dα2y=0. If y=0, the tumor-free equilibrium is at ℰ¯0=(x¯,y¯)=(σ/θ,0). This steady state is always exist, since σ/θ>0. From the analysis, it is easy to prove that the tumor-free equilibrium ℰ¯0=(σ/θ,0) of the model (10) is asymptotically stable if threshold parameter (the minimum tumor-clearance parameter) ℛ¯0=aθ/σ<1 and unstable if ℛ¯0>1.However, ify≠0, the steady states are obtained by solving
(11)ωaby2-a(ω+θb)y+aθ-σ=0.
In this case, we have two endemic equilibria, ℰ¯1 and ℰ¯2:
(12)ℰ¯1=(x¯1,y¯1),wherex¯1=-a(bθ-ω)-Δ2ω,ℰ¯1=(x¯1,y¯1),wherex¯1y¯1=a(bθ+ω)+Δ2abω,ℰ¯2=(x¯2,y¯2),wherex¯2=-a(bθ-ω)+Δ2ω,ℰ¯2=(x¯2,y¯2),wherex¯2y¯2=a(bθ+ω)-Δ2abω
with Δ=a2(bθ-ω)2+4σωab>0. The Jacobian matrix of the system (10) at the endemic equilibrium ℰ¯1 is
(13)J(ℰ¯1)=(ωy¯1-θωx¯1-y¯1a-2aby¯1-x¯1).Proposition 4.
Assume that the endemic equilibriumℰ¯1 exists and has nonnegative coordinates. If ℛ¯0=aθ/σ<1, then tr(J(ℰ¯1))>0 and ℰ¯1 is unstable.Proof.
Since(14)tr(J(ℰ¯1))=ω2-ω(ab+bθ)-ab2θ2bωtr(J(ℰ¯1))=+ω-ab2abωa2(bθ+ω)2-4abω(aθ-σ),
then inequality tr(J(ℰ¯1))>0 is true if
(15)a[ω2-ω(ab+bθ)-ab2θ]>(ab-ω)a2(bθ+ω)2-4abω(aθ-σ).
Therefore, when aθ<σ, we have ω2-ωb(a+θ)-aθb2>0, and hence both sides of the inequality are positive. Therefore if the equilibrium point ℰ¯1 exists and has nonnegative coordinates, then tr(J(ℰ¯1))>0 and the point (ℰ¯1) is unstable, whenever ℛ¯0=aθ/σ<1.Similarly, we arrive at the following proposition.Proposition 5.
If the pointℰ¯2 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:
(16)Dα1H=s-μHH+rH(1-H+IHmax)-k1VH,Dα2I=k1′VH-μII,Dα3V=MμbI-k1VH-μVV,0.5<αi≤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, μH 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. k1′ is the rate at which infected cells become actively infected. μI 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, μb 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, μV is the loss rate of virus. The initial conditions for infection by free virus are H(0)=H0, I(0)=I0, and V(0)=V0.Theorem 6.
According to Remark3 and the fact that s≥0, then the system (16) has a unique solution (H,I,V)T which remains in ℝ+3 and bounded by Hmax; [11].


### 3.1. Equilibria and Local Stability of Model (16)

To evaluate the equilibrium points of system (16), we put Dα1H(t)=Dα2I(t)=Dα3V(t)=0. Then the infection-free equilibrium point (uninfected steady state) is ℰ0*=(H0,0,0), and endemic equilibrium point (infected steady state) is ℰ+*=(H*,I*,V*), where
(17)H*=μVμIk1′Mμb-k1μI,I*=k1′H*V*μI,V*=μI[(s+(r-μH)H*)Hmax-rH*2]H*[k1′rH*-k1μIHmax].
The Jacobian matrix J(ℰ0*) for system (16) evaluated at the uninfected steady state ℰ0* is then given by
(18)J(ℰ0*)=(-μH+r-2rH0Hmax-rH0Hmax-k1H00-μIk1′H00Mμb-(k1H0+μV)).
Let us introduce the following definition and assumption to ease the analysis.Definition 7.
The threshold parameterℛ0* (the minimum infection-free parameter) is the parameter that has the property that if ℛ0*<1, then the endemic infected state does not exist, while if ℛ0*>1, the endemic infected state persists, where
(19)ℛ0*=k1′μbMH0μI(μV+k1H0).We also assume that(20)H0={r-μH+[(r-μH)2+4rsHmax-1]1/2}2rHmax-1.The uninfected steady state is asymptotically stable if all of the eigenvaluesλ of the Jacobian matrix J(ℰ0*), given by (18), have negative real parts. The characteristic equation det(J(ℰ0)-I)=0 becomes
(21)(λ+μH-r+2rH0Hmax)(λ2+Bλ+C)=0,
where B=μI+kIH0+μV and C=μI(k1H0+μV)-k1′μbMH0. Hence, the three roots of the characteristic equation (21) are
(22)λ1=-μH+r-rH0Hmax≡-(r-μH)2+4rsHmax-1<0,λ2,3=12[-B±B2-4C].Proposition 8.
Ifℛ0*≡(k1′μbMH0)/μI(μV+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ℰ0*, one has three cases. (i)
Ifℛ0*<0, the uninfected state is asymptotically stable and the infected steady state ℰ+ does not exist (unphysical).(ii)
Ifℛ0*=1, then C=0, and from (21) implies that one eigenvalue λ=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ℛ0*>1, then C<0, and thus at least one eigenvalue will be positive real root. Thus, the uninfected state ℰ0 is unstable, and the endemically infected state ℰ+* emerges.To study the local stability of the positive infected steady statesℰ+* for ℛ0*>1, we consider the linearized system of (16) at ℰ+*. The Jacobian matrix at ℰ+* becomes
(23)J(ℰ+*)=(-L*-rH*Hmax-k1H*k1′V*-μIk1′H*k1V*Mμb-(k1H*+μV)).
Here
(24)L*=-[μH-r+k1V*+r(2H*+I*)Hmax].
Then the characteristic equation of the linearized system is
(25)P(λ)=λ3+a1λ2+a2λ+a3=0,(26)a1=μI+μV+k1H*+L*,a2=L*(μI+μV+k1H*)+μI(μV+k1H*)-k12H*V*-k1′H*(Mμb-rV*Hmax),a3=k1′H*[k1MμbV*-rμVV*Hmax-L*Mμb]+L*μI(μV+k1H*)-μIk12H*V*.
The infected steady state ℰ+* 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.


## 3.1. Equilibria and Local Stability of Model (16)

To evaluate the equilibrium points of system (16), we put Dα1H(t)=Dα2I(t)=Dα3V(t)=0. Then the infection-free equilibrium point (uninfected steady state) is ℰ0*=(H0,0,0), and endemic equilibrium point (infected steady state) is ℰ+*=(H*,I*,V*), where
(17)H*=μVμIk1′Mμb-k1μI,I*=k1′H*V*μI,V*=μI[(s+(r-μH)H*)Hmax-rH*2]H*[k1′rH*-k1μIHmax].
The Jacobian matrix J(ℰ0*) for system (16) evaluated at the uninfected steady state ℰ0* is then given by
(18)J(ℰ0*)=(-μH+r-2rH0Hmax-rH0Hmax-k1H00-μIk1′H00Mμb-(k1H0+μV)).
Let us introduce the following definition and assumption to ease the analysis.Definition 7.
The threshold parameterℛ0* (the minimum infection-free parameter) is the parameter that has the property that if ℛ0*<1, then the endemic infected state does not exist, while if ℛ0*>1, the endemic infected state persists, where
(19)ℛ0*=k1′μbMH0μI(μV+k1H0).We also assume that(20)H0={r-μH+[(r-μH)2+4rsHmax-1]1/2}2rHmax-1.The uninfected steady state is asymptotically stable if all of the eigenvaluesλ of the Jacobian matrix J(ℰ0*), given by (18), have negative real parts. The characteristic equation det(J(ℰ0)-I)=0 becomes
(21)(λ+μH-r+2rH0Hmax)(λ2+Bλ+C)=0,
where B=μI+kIH0+μV and C=μI(k1H0+μV)-k1′μbMH0. Hence, the three roots of the characteristic equation (21) are
(22)λ1=-μH+r-rH0Hmax≡-(r-μH)2+4rsHmax-1<0,λ2,3=12[-B±B2-4C].Proposition 8.
Ifℛ0*≡(k1′μbMH0)/μI(μV+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ℰ0*, one has three cases. (i)
Ifℛ0*<0, the uninfected state is asymptotically stable and the infected steady state ℰ+ does not exist (unphysical).(ii)
Ifℛ0*=1, then C=0, and from (21) implies that one eigenvalue λ=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ℛ0*>1, then C<0, and thus at least one eigenvalue will be positive real root. Thus, the uninfected state ℰ0 is unstable, and the endemically infected state ℰ+* emerges.To study the local stability of the positive infected steady statesℰ+* for ℛ0*>1, we consider the linearized system of (16) at ℰ+*. The Jacobian matrix at ℰ+* becomes
(23)J(ℰ+*)=(-L*-rH*Hmax-k1H*k1′V*-μIk1′H*k1V*Mμb-(k1H*+μV)).
Here
(24)L*=-[μH-r+k1V*+r(2H*+I*)Hmax].
Then the characteristic equation of the linearized system is
(25)P(λ)=λ3+a1λ2+a2λ+a3=0,(26)a1=μI+μV+k1H*+L*,a2=L*(μI+μV+k1H*)+μI(μV+k1H*)-k12H*V*-k1′H*(Mμb-rV*Hmax),a3=k1′H*[k1MμbV*-rμVV*Hmax-L*Mμb]+L*μI(μV+k1H*)-μIk12H*V*.
The infected steady state ℰ+* 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(27)DαX(t)=F(t,X(t),P),t∈[0,T],0<α≤1,X(0)=X0.
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
(28)∥F(t,X(t),P)-F(t,Y(t),P)∥≤K∥X(t)-Y(t)∥,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaK>0,
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 KTα/Γ(α+1)<1.Proof.
Using the definitions of Section1, 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
(29)X(t)=X(0)+1Γ(α)∫0t‍(t-s)α-1F(s,X(s),P)ds
which also is a Volterra equation of the second kind. Define the operator ℒ, such that
(30)ℒX(t)=X(0)+1Γ(α)∫0t‍(t-s)α-1F(s,X(s),P)ds.
Then, we have
(31)∥ℒX(t)-ℒY(t)∥≤1Γ(α)∫0t‍(t-s)α-1∥F(s,X(s),P)-F(s,Y(s),P)∥ds≤KΓ(α)∫0t‍(t-s)α-1sups∈[0,T]|X(s)-Y(s)|ds≤KΓ(α)∥X-Y∥∫0t‍sα-1ds≤KΓ(α)∥X-Y∥Tα.
Thus for KTα/Γ(α+1)<1, we have
(32)∥ℒX(t)-ℒY(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 𝒯={t0,t1,…,tN}, such that t0=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 α, 0<α≤1, and using implicit Euler’s approximation as follows (see [44]):
(33)D*αxi(tn)=1Γ(1-α)∫0t‍dxi(s)ds(tn-s)-αds≈1Γ(1-α)∑j=1n‍∫(j-1)hjh‍[xij-xij-1h+O(h)](nh-s)-αds=1(1-α)Γ(1-α)×∑j=1n{[xij-xij-1h+O(h)]‍×[(n-j+1)1-α-(n-j)1-α][xij-xij-1h+O(h)]}h1-α=1(1-α)Γ(1-α)1hα×∑j=1n‍[xij-xij-1][(n-j+1)1-α-(n-j)1-α]+1(1-α)Γ(1-α)×∑j=1n‍[xij-xij-1][(n-j+1)1-α-(n-j)1-α]O(h2-α).
Setting
(34)𝒢(α,h)=1(1-α)Γ(1-α)1hα,ωjα=j1-α-(j-1)1-α,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa(whereω1α=1),
then the first-order approximation method for the computation of Caputo’s fractional derivative is then given by the expression
(35)D*αxi(tn)=𝒢(α,h)∑j=1n‍ωjα(xin-j+1-xin-j)+O(h).
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(36)D*αu(t)=ρ0u(t)+ρ1,U(0)=U0,
such that 0<α≤1, and ρ0<0, ρ1>0 are constants.Theorem 12.
The fully implicit numerical approximation (35), 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
(37)(1-ρ0𝒢α,h)ζn=ζn-1+∑j=2n‍ωj(α)(ζn-j-ζn-j+1)+ρ1𝒢α,h,n≥2
or
(38)ζn=ζn-1+∑j=2nωj(α)(ζn-j-ζn-j+1)+ρ1/𝒢α,h(1-(ρ0/𝒢α,h)),n≥2.
Since (1-(ρ0/𝒢α,h))≥1 for all 𝒢α,h, then
(39)ζ1≤ζ0,(40)ζn≤ζn-1+∑j=2n‍ωj(α)(ζn-j-ζn-j+1),n≥2.
Thus, for n=2, the previous inequality implies
(41)ζ2≤ζ1+ω2(α)(ζ0-ζ1).
Using the relation (39) and the positivity of the coefficients ω2, we get
(42)ζ2≤ζ1.
Repeating the process, we have from (40)
(43)ζn≤ζn-1+∑j=2n‍ωj(α)(ζn-j-ζn-j+1)≤ζn-1,
since each term in the summation is negative. Thus ζn≤ζn-1≤ζn-2≤⋯≤ζ0. With the assumption that ζn=|Un|≤ζ0=|U0|, which entails ∥Un∥≤∥U0∥, 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<α≤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 α=1) and fractional-order (α<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.Figure 1
Numerical simulations of the FODEs model (5) when α=0.75 and α=0.95 with (a=r1=r2=1; d1=0.3; d2=0.7; k1=0.3; k2=0.7) and when d1=0.7; d2=0.3. The system converges to the stable steady states ℰ1, ℰ2. The fractional derivative damps the oscillation behavior.Figure 2
Numerical simulations of the FODEs model (10) when s=0.1181, ω=0.1184, θ=0.1747, r=0.636, and b=0.002. α1=α2=1,0.9,0.7. The endemic state ℰ¯2 is locally asymptotically stable when ℛ¯0=aθ/σ>1. The fractional derivative damps the oscillation behavior.Figure 3
Numerical simulations of the FODEs model (10) when s=0.1181, ω=0.1184, θ=0.3747, r=1.636, and b=0.002. α1=α2=1,0.9,0.7. The infection-free steady state ℰ¯0 is locally asymptotically stable when ℛ¯0=aθ/σ<1.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 α, 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, the solution of the fractional order model tends to the fixed point over a longer period of time.Figure 4
Numerical simulations of the FODEs model (16) when s=20day-1mm-3, μH=0.02day-1, μI=0.26day-1, μb=0.24day-1, μV=2.4day-1, k1=2.4×10-5mm3day-1, k1′=2×10-5mm3day-1, r=0.3day-1, Hmax=1500mm-3, and M=1400. The trajectories of the system approach to the steady state ℰ+*.


## 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(36)D*αu(t)=ρ0u(t)+ρ1,U(0)=U0,
such that 0<α≤1, and ρ0<0, ρ1>0 are constants.Theorem 12.
The fully implicit numerical approximation (35), 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
(37)(1-ρ0𝒢α,h)ζn=ζn-1+∑j=2n‍ωj(α)(ζn-j-ζn-j+1)+ρ1𝒢α,h,n≥2
or
(38)ζn=ζn-1+∑j=2nωj(α)(ζn-j-ζn-j+1)+ρ1/𝒢α,h(1-(ρ0/𝒢α,h)),n≥2.
Since (1-(ρ0/𝒢α,h))≥1 for all 𝒢α,h, then
(39)ζ1≤ζ0,(40)ζn≤ζn-1+∑j=2n‍ωj(α)(ζn-j-ζn-j+1),n≥2.
Thus, for n=2, the previous inequality implies
(41)ζ2≤ζ1+ω2(α)(ζ0-ζ1).
Using the relation (39) and the positivity of the coefficients ω2, we get
(42)ζ2≤ζ1.
Repeating the process, we have from (40)
(43)ζn≤ζn-1+∑j=2n‍ωj(α)(ζn-j-ζn-j+1)≤ζn-1,
since each term in the summation is negative. Thus ζn≤ζn-1≤ζn-2≤⋯≤ζ0. With the assumption that ζn=|Un|≤ζ0=|U0|, which entails ∥Un∥≤∥U0∥, 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<α≤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 α=1) and fractional-order (α<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.Figure 1
Numerical simulations of the FODEs model (5) when α=0.75 and α=0.95 with (a=r1=r2=1; d1=0.3; d2=0.7; k1=0.3; k2=0.7) and when d1=0.7; d2=0.3. The system converges to the stable steady states ℰ1, ℰ2. The fractional derivative damps the oscillation behavior.Figure 2
Numerical simulations of the FODEs model (10) when s=0.1181, ω=0.1184, θ=0.1747, r=0.636, and b=0.002. α1=α2=1,0.9,0.7. The endemic state ℰ¯2 is locally asymptotically stable when ℛ¯0=aθ/σ>1. The fractional derivative damps the oscillation behavior.Figure 3
Numerical simulations of the FODEs model (10) when s=0.1181, ω=0.1184, θ=0.3747, r=1.636, and b=0.002. α1=α2=1,0.9,0.7. The infection-free steady state ℰ¯0 is locally asymptotically stable when ℛ¯0=aθ/σ<1.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 α, 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, the solution of the fractional order model tends to the fixed point over a longer period of time.Figure 4
Numerical simulations of the FODEs model (16) when s=20day-1mm-3, μH=0.02day-1, μI=0.26day-1, μb=0.24day-1, μV=2.4day-1, k1=2.4×10-5mm3day-1, k1′=2×10-5mm3day-1, r=0.3day-1, Hmax=1500mm-3, and M=1400. The trajectories of the system approach to the steady state ℰ+*.


## 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 ℛ0 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α 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 the previous states. The analysis can be extended to other models related to the immune response systems.



---

*Source: 816803-2013-08-07.xml*