Numerical Modeling of Fractional-Order Biological Systems

and Applied Analysis 3 the nonbiological interior solution where both immune effectors coexist. The equilibrium points of the system (5) are E 0 = (0, 0, 0) , E 1 = (√ d 1 k 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][2][3][4][5], physics [6][7][8][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][19][20][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.

Abstract and Applied Analysis
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 1 = 1 [ , ] be the class of Lebesgue integrable functions on [ , ], < < ∞.
Definition 1. The fractional integral (or the Riemann-Liouville integral) of order ∈ R + of the function ( ), > 0 ( : The fractional derivative of order ∈ ( − 1, ) of ( ) is defined by two (nonequivalent) ways: (i) Riemann-Liouville fractional derivative: take fractional integral of order ( − ), and then take th derivative as follows: (ii) Caputo-fractional derivative: take th derivative, and then take a fractional integral of order ( − ) We notice that the definition of time-fractional derivative of a function ( ) at = involves an integration and calculating time-fractional derivative that requires all the past history, that is, all the values of ( ) from = 0 to = . 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][30][31]).

Remark 2.
Let , ∈ R + and ∈ (0, 1). Then The generalized mean value theorem and another property are defined in the following remark [32].

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][34][35][36][37], but fractionalorder 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 cancerimmune system interaction. The model includes two immune effectors: 1 ( ), 2 ( ) (such as cytotoxic T cells and natural killer cells), interacting with the cancer cells, ( ), 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 where ≡ ( ), 1 ≡ 1 ( ), 2 ≡ 2 ( ), and , 1 , 2 , 1 , 2 , 1 , and 2 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 ( 1 1 /(1 − 1 )) ≪ ( 2 2 /(1 − 2 )) to avoid the nonbiological interior solution where both immune effectors coexist. The equilibrium points of the system (5) are The first equilibrium E 0 is the nave, the second E 1 is the memory, and the third E 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 1 < 2 and 1 < 1. While the endemic state is locally asymptotically stable if 2 < 1 and 2 < 1, there is bifurcation at 1 = 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, ( ); the tumor cells, ( ); and the concentration of IL-2 in the single tumor-site compartment, ( ). We consider the classic bilinear model that includes Holling function of type I and external effector cells 1 and input of IL-2, 2 . (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: with initial conditions: (0) = 0 , (0) = 0 , and (0) = 0 . The parameter 1 is the antigenicity rate of the tumor (immune response to the appearance of the tumor), and 1 is the external source of the effector cells, with rate of death 2 . The parameter 4 incorporates both multiplication and death of tumor cells. The maximal carrying capacity of the biological environment for tumor cell is −1 5 ; 3 is considered as the cooperation rate of effector cells with interleukin-2 parameter; 6 is the rate of tumor cells; 7 is the competition rate between the effector cells and the tumor cells. External input of IL-2 into the system is 2 , and the rate loss parameter of IL-2 cells is 8 .
In the absence of immunotherapy with IL-2, we have 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 Therefore, after the previous substitution into (8) and replacing by , the model becomes (10). The steady states of the reduced model (10) are again the intersection of the null clines 1 = 0, 2 = 0. If = 0, the tumorfree equilibrium is at E 0 = ( , ) = ( / , 0). This steady state is always exist, since / > 0. From the analysis, it is easy to prove that the tumor-free equilibrium E 0 = ( / , 0) of the model (10) is asymptotically stable if threshold parameter (the minimum tumor-clearance parameter) R 0 = / < 1 and unstable if R 0 > 1. However, if ̸ = 0, the steady states are obtained by solving

Equilibria and Local Stability of Model
In this case, we have two endemic equilibria, E 1 and E 2 : with Δ = 2 ( − ) 2 + 4 > 0. The Jacobian matrix of the system (10) at the endemic equilibrium E 1 is Proposition 4. Assume that the endemic equilibrium E 1 exists and has nonnegative coordinates. If R 0 = / < 1, then tr( (E 1 )) > 0 and E 1 is unstable.
then inequality tr( (E 1 )) > 0 is true if Therefore, when < , we have 2 − ( + )− 2 > 0, and hence both sides of the inequality are positive. Therefore if the equilibrium point E 1 exists and has nonnegative coordinates, then tr( (E 1 )) > 0 and the point ( Similarly, we arrive at the following proposition.

Proposition 5. If the point E 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.

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: Here = ( ) represents the concentration of healthy CD4 + T cells at time , = ( ) represents the concentration of infected CD4 + T cells, and = ( ) is the concentration of free HIV at time . In this system, the logistic growth of the healthy CD4 + T-cells is (1 − (( + )/ max )), and the proliferation of infected cells is neglected. The parameter is the source of CD4 + T cells from precursors, is the natural death rate of CD4 + T cells ( max > , cf. [39, page 85]), is their growth rate (thus, > in general), and max is their carrying capacity. The parameter 1 represents the rate of infection of T cells with free virus. 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 viral particles are released by each lysing cell, this term is multiplied by the parameter 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 (0) = 0 , (0) = 0 , and (0) = 0 . Theorem 6. According to Remark 3 and the fact that ≥ 0, then the system (16) has a unique solution ( , , ) which remains in R 3 + and bounded by max ; [11]. (16). To evaluate the equilibrium points of system (16)

Equilibria and Local Stability of Model
The Jacobian matrix (E * 0 ) for system (16) evaluated at the uninfected steady state E * 0 is then given by ) .
Let us introduce the following definition and assumption to ease the analysis.
Definition 7. The threshold parameter R * 0 (the minimum infection-free parameter) is the parameter that has the property that if R * 0 < 1, then the endemic infected state does not exist, while if R * 0 > 1, the endemic infected state persists, where Abstract and Applied Analysis 5 We also assume that The uninfected steady state is asymptotically stable if all of the eigenvalues of the Jacobian matrix (E * 0 ), given by (18), have negative real parts. The characteristic equation det( (E 0 ) − ) = 0 becomes where = + 0 + and = ( 1 0 + )− 1 0 . Hence, the three roots of the characteristic equation (21) are Proposition 8. If R * 0 ≡ ( 1 0 )/ ( + 1 0 ) < 1, then > 0 and the three roots of the characteristic equation (21) will have negative real parts.

Corollary 9.
In case of uninfected steady state E * 0 , one has three cases.
(i) If R * 0 < 0, the uninfected state is asymptotically stable and the infected steady state E + does not exist (unphysical).
(ii) If R * 0 = 1, then = 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 R * 0 > 1, then < 0, and thus at least one eigenvalue will be positive real root. Thus, the uninfected state E 0 is unstable, and the endemically infected state E * + emerges.
To study the local stability of the positive infected steady states E * + for R * 0 > 1, we consider the linearized system of (16) at E * + . The Jacobian matrix at E * + becomes ) . Here Then the characteristic equation of the linearized system is The infected steady state E * + 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, 1 > 0, 3 > 0, and 1 2 > 3 .

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 where ( ) is the solution of the perturbed system. 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 which also is a Volterra equation of the second kind. Define the operator L, such that Then, we have Thus for /Γ( + 1) < 1, we have 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 = { 0 , 1 , . . . , }, such that 0 = 0 and = . 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]): * ( ) Setting then the first-order approximation method for the computation of Caputo's fractional derivative is then given by the 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.

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 ℎ → 0.
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].

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 ℎ = 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 fractionalorder ( < 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 , 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, Abstract and Applied Analysis

Conclusions
In fact, fractional-order differential equations are generalizations of integer-order differential equations. Using fractionalorder 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 R 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 * is continuously depends on all the previous states at ≤ * . 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 diseasefree 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.