Abstract and Applied Analysis

Volume 2013 (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

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

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

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.

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 (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.

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 be the class of Lebesgue integrable functions on , .

Definition 1. The fractional integral (or the Riemann-Liouville integral) of order of the function , () is defined by

The fractional derivative of order 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 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–31]).

Remark 2. Let and . Then

(i)if and , then ; (ii) uniformly on , , where ; (iii) weakly; (iv)if is absolutely continuous on , then ; (v)thus (Riemann-Liouville sense) and (Caputo sense).

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

Remark 3.

(i)Suppose that and for , and then we have

(ii)If (i) holds and for all , then is nondecreasing for each . If for all , then is nonincreasing for each .

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: , (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 , , , and , , , , , , and 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 to avoid the nonbiological interior solution where both immune effectors coexist. The equilibrium points of the system (5) are

The first equilibrium is the nave, the second is the memory, and the third 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 and . While the endemic state is locally asymptotically stable if and , there is bifurcation at . 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 and input of IL-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: , , and . The parameter 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 . The parameter incorporates both multiplication and death of tumor cells. The maximal carrying capacity of the biological environment for tumor cell is ; is considered as the cooperation rate of effector cells with interleukin-2 parameter; is the rate of tumor cells; is the competition rate between the effector cells and the tumor cells. External input of IL-2 into the system is , and the rate loss parameter of IL-2 cells is .

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

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 , . If , the tumor-free equilibrium is at . This steady state is always exist, since . From the analysis, it is easy to prove that the tumor-free equilibrium of the model (10) is asymptotically stable if threshold parameter (the minimum tumor-clearance parameter) and unstable if .

However, if , the steady states are obtained by solving

In this case, we have two endemic equilibria, and :

with . The Jacobian matrix of the system (10) at the endemic equilibrium is

Proposition 4. Assume that the endemic equilibrium exists and has nonnegative coordinates. If , then and is unstable.

Proof. Since

then inequality is true if

Therefore, when , we have , and hence both sides of the inequality are positive. Therefore if the equilibrium point exists and has nonnegative coordinates, then and the point is unstable, whenever .

Similarly, we arrive at the following proposition.

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

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 , 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 (, cf. [39, page 85]), is their growth rate (thus, in general), and is their carrying capacity. The parameter represents the rate of infection of T cells with free virus. 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 , , and .

Theorem 6. According to Remark 3 and the fact that , then the system (16) has a unique solution which remains in and bounded by ; [11].

3.1. Equilibria and Local Stability of Model (16)

To evaluate the equilibrium points of system (16), we put . Then the infection-free equilibrium point (uninfected steady state) is , and endemic equilibrium point (infected steady state) is , where

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

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

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

We also assume that

The uninfected steady state is asymptotically stable if all of the eigenvalues of the Jacobian matrix , given by (18), have negative real parts. The characteristic equation becomes

where and . Hence, the three roots of the characteristic equation (21) are

Proposition 8. If , then and the three roots of the characteristic equation (21) will have negative real parts.

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

(i)If , the uninfected state is asymptotically stable and the infected steady state does not exist (unphysical). (ii)If , then , and from (21) implies that one eigenvalue and the remaining two eigenvalues have negative real parts. The uninfected and infected steady states collide, and there is a transcritical bifurcation. (iii)If , then , and thus at least one eigenvalue will be positive real root. Thus, the uninfected state is unstable, and the endemically infected state emerges.

To study the local stability of the positive infected steady states for , we consider the linearized system of (16) at . The Jacobian matrix at becomes

Here

Then the characteristic equation of the linearized system is

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, , , and .

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

Here , is the set of parameters appear in the model, and satisfies the Lipschitz condition

where is the solution of the perturbed system.

Theorem 10. Problem (27) has a unique solution provided that the Lipschitz condition (28) is satisfied and .

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 , such that

Then, we have

Thus for , 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 , such that 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 , , 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 expression

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 .

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

such that , and , are constants.

Theorem 12. The fully implicit numerical approximation (35), to test problem (36) for all , is consistent and unconditionally stable.

Proof. We assume that the approximate solution of (36) is of the form , and then (36) can be reduced to

or

Since for all , then

Thus, for , the previous inequality implies

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

Repeating the process, we have from (40)

since each term in the summation is negative. Thus . With the assumption that , which entails , 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 and 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 ) and fractional-order () 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 and with (; ; ; ; ) and when ; . The system converges to the stable steady states , . The fractional derivative damps the oscillation behavior.

Figure 2: Numerical simulations of the FODEs model (10) when , , , , and . . The endemic state is locally asymptotically stable when . The fractional derivative damps the oscillation behavior.

Figure 3: Numerical simulations of the FODEs model (10) when , , , , and . . The infection-free steady state is locally asymptotically stable when .

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 , , , , , , , , , and . 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 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 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.

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

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.

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.

K. S. Cole, “Electric conductance of biological systems,” Cold Spring Harbor Symposia on Quantitative Biology, pp. 107–116, 1993.

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.

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.

L. Debnath, “Recent applications of fractional calculus to science and engineering,” International Journal of Mathematics and Mathematical Sciences, no. 54, pp. 3413–3442, 2003.

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.

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

G. M. Zaslavsky, “Chaos, fractional kinetics, and anomalous transport,” Physics Reports, vol. 371, no. 6, pp. 461–580, 2002.

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.

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.

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

K. Assaleh and W. M. Ahmad, “Modeling of speech signals using fractional calculus,” in Proceedings of the 9th International Symposium on Signal Processing and its Applications (ISSPA '07), Sharjah, United Arab Emirates, February 2007.

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.

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.

A. Rocco and B. J. West, “Fractional calculus and the evolution of fractal phenomena,” Physica A, vol. 265, no. 3, pp. 535–546, 1999.

V. D. Djordjević, J. Jarić, B. Fabry, J. J. Fredberg, and D. Stamenović, “Fractional derivatives embody essential features of cell rheological behavior,” Annals of Biomedical Engineering, vol. 31, no. 6, pp. 692–699, 2003.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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

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

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

C. Li and Y. Ma, “Fractional dynamical system and its linearization theorem,” Nonlinear Dynamics, vol. 71, no. 4, pp. 621–633, 2013.

Z. M. Odibat and N. T. Shawagfeh, “Generalized Taylor's formula,” Applied Mathematics and Computation, vol. 186, no. 1, pp. 286–293, 2007.

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.

A. Gökdoğan, 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.

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.

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.

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.

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.

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.

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. 183–199, 2000.

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.

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.

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.

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.

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.