dynamical analysis of mathematical model presented by fractional differential equations, describing the interaction between leukemic cancer cells, t cells and drug treatment with a drug optimal control

Research Article

黑料网

G H Erjaee^1*, M Shahbazi¹ and A Erjaee²

¹Mathematics Department, Shiraz University, Shiraz, Iran

²Assistant Professor, Pediatric Department, Shiraz University, Medical School, Shiraz, Iran

^*Corresponding authors:

G. H. Erjaee
Mathematics Department
Shiraz University
Shiraz, Iran
E-mail: erjaee@shirazu.ac.ir

Received August 21, 2011; Published October 01, 2012

Citation: Erjaee GH, Shahbazi M, Erjaee A (2012) Dynamical Analysis of Mathematical Model Presented by Fractional Differential Equations, Describing the Interaction Between Leukemic Cancer Cells, T Cells and Drug Treatment with a Drug Optimal Control. 1:390. doi:10.4172/scientificreports.390

Copyright: © 2012 Erjaee GH, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Obviously, one's primary motivation in producing any mathematical model is to describe a natural or artificial phenomenon by means of a model equation whose behavior is as close as possible to that of the original phenomenon. But this is often difficult, particularly when we are dealing with nonlinear behavior in natural complex phenomena such as the interaction between variety of cancers and immune cells. Thus, the choice of modeling schemes that produce models whose dynamics resembles those of their physical counterparts is a major challenge in mathematical molding. Ordinary differential equations and partial differential equations are just some of the mathematical tools that are been used in deriving these mathematical models. In this article, we have used fractional differential equations for our derivation. In this derivation, a mathematical model describes the growth or terminates myelogenous leukemia blood cancer’s cells against naive T-cell and effective T-cell cells of body. Using this model, we have studied the dynamic behavior describing the transaction between bodies' effective T cell, naive T cell and chronic myelogenous leukemia in one side and drug in other side. The most important feature of equations with fractional order derivatives is their non-localization. We expect that our fractional differential equations model will be superior to its ordinary differential equations counterpart in facilitating understanding of the natural immune interactions to tumor and of the detrimental side-effects which chemotherapy may have on a patient’s immune system. Using this system, we will study the optimized drug dose in chronic myelogenous leukemia treatment with two methods namely targeted therapy and broad cytotoxic therapy.

Keywords

Non-linear dynamics; Fractional differential equations; Cancer chemotherapy; Optimal drug control

Mathematical Subject Classification

37N25, 37N35 and 44A55.

Introduction

Drug therapy is a way to cure cancer. It is clear that, the drug dose is important for cancer specialists. Considering the weakness of immunology system in cancer affected patients, drug overdose may results in additional problems for their body.

Chronic myelogenous leukemia (CML) is a kind of blood cancer and occurs in adults about 15 percent [1]. The age average for blood cancer patients ranges between 45-55 years old. The occurrence rate is one or three among per 100000 individuals [1]. In most CML types, leukemic cells are divided to an unmoral Chromosome which can't be found neither in non-blood white cells nor in any of body’s cells [2]. Consequently, a displacement will be occurred between Chromosome 9 and 22 in such a way that Chromosome 9 will be longer than normal state and Chromosome 22 will be short which is called Chromosome Philadelphia (Ph) [2]. This Chromosome produce an unmoral protein namely Tyrosine kinase (Ber-Abl). This change results in conversion of bone marrow cells to abnormal leukemic cells [3,4]. Structural knowledge of Ber-Abl has led to the development of a drug, imatinib mesylate (known as Gleevec in the U.S. and henceforth referred to as imatinib), that blocks the abnormal protein, thus removing the proliferative advantage that it provides to cancer cells. Imatinib is highly specific to binding with Ber-Abl, and hence, generally has mild side effects [2]. One of the obtained successes in this regard was related to Food and Drug Administration (FDA) for detecting CML diseases in December, 2002 [5]. Bone marrow transplantation is a curative treatment option for some patients, but transplant-related mortality rates can be above 40 percent [6]. Prior to the development of imatinib, treatments such as hydroxyuea, cytarabine, interferonalfa or combination of them were used to treat CML [5]. The action of these therapies is against broad classes of cells, and so treatment usually results in severe side effects. Several studies suggest that combination of imatinib with a broader chemotherapy has the potential to perform better than imatinib alone [7-10].

Recently, the models used for analyzing the cancer reaction against drug therapy could assist physicians in cancer treatment. Therefore, using optimized control methods which minimized damages to body can play an important role in cancer treatment. In this field, researchers such as Fokas et al. [11] and Adimy et al. [12] have presented CML models in 1991 and 2005, respectively. A mathematical model has also presented by Afenya and Bentil in 1998 [13] for blood cancer. Periodic mathematical models for CML have presented by Mackey and Menjouet in 2004.

In this article, the ordinary differential equation (ODE) which is presented by Moore and Li for brain blood cancer [14] is re-derived using fractional differential equation (FDE). We claim that our FDE model will be superior to its ODE counterpart in facilitating understanding of the natural immune interactions to tumor and of the detrimental side effects which chemotherapy may have on a patient’s immune system. Summarizing the advantages of our FDE model over previous ones, first as Andrew Einstein described in his research [15] at Mount Sinai School of Medicine, some cells in various body organs have a rugged surface that cannot be properly understand using ordinary calculus, but it may be amenable to studies using fractional calculus. In addition, since differentiability is not required in accordance with the definition of FDE with order between 0 and 1, so these equations can be used in non-smooth domains. Furthermore, whereas in the definition of the classical derivative of a function at a point we use just two points in the neighborhood of that point in the definition of the fractional derivative we use all the points in a neighborhood of the point. Obviously, by using all of the information available in these points we obtain more accurate results in subsequent applications. The last property that is so-called non-local property will closely reflect reality and is a primary reason why FDE are increasingly applied to dynamical systems.

Therefore, in this article, first we will introduce a FDE model to present the interaction between naive T cells, effectors T cells, and CML cancer cells in cancer dormancy. Then, we will discuss the dynamical behavior of this model by identifying the fixed points and determining their stability characteristics. To find the solutions of this FDE system, we will discretize the system by using Grunwald- Letnikov discretization method [16,17] then, we will find the results by using software tools such as MATLAB™. In this FDE model, by adding chemotherapy drug concentration to the interaction between naive T cells, effectors T cells, and CML cancer cells and considering the same three cells populations as in the first FDE, we will have our second model in the form of FDE. Now, similar to the way in which have been done in classical ODE systems, we will discuss the dynamic behavior of the first system and determine the stability type of the various feasible fixed points. One of the main goals in using fractional order instead of classical integer order derivative in our model is to obtain more accurate results in chemotherapy optimal control. For this optimality, similar to the targeted therapy (such as imatinib) and broad cytotoxic therapy (such as cytarabine) methods used by Moore and Li [14], we will use the processors in our FDE model. Obviously, in the processing of this optimality we need to solve our FDE system numerically. To facilitate this solution, as in above first FDE system, we will apply Grunwald-Letnikov method to discretize the model and then, use MATLAB™ software to find the results.

As we stated before, we will expect more accurate results in solving our FDE systems as compared to the results found by classical ODE methods.

Dynamical Analyses of Tumor without Treatment in the FDE Model

The first model that we consider here is a three cells population model describing the interaction between the cancer cell population (C), the naive T cell population (T_n) and effector T cell population (T_e). We assume that the effector T cells are specific to CML, activated by the presence of CML antigen. If we suppose these three cells evolve with independent variable time, then we can present our model in the form of FDE as follows.

This model is similar to the classical ODE model presented by Seema Nanda and Helen Moore [14]. Here, we use the same order derivatives α∈(0,1) for all three equations, where T_n(0), T_e (0) and C(0) are known initial values. All of the parameter values in the above equations are assumed to be positive. The structure of the equations guarantees non-negative solutions for the state variables, (t), (t) and C(t) [2]. The negative terms in the above equations represent losses from the cell populations while the positive terms are source terms for the cell populations (Figure 1). The last term in the first equation (which has a Michaelis–Menten factor) represents the decline in T_n cells due to encounters with CML antigen in the lymph nodes. As the population of T_n cells is very small in comparison to the CML cells, this term takes into account the saturation effect of CML cells in the lymph nodes. Since some of these lost T_n cells reappear as effector T cells (T_e), a Michaelis-Menten factor also shows up in the first term in Eq. (2). The model assumes Gompertz growth for CML cells as indicated by the growth term in Eq. (3). Also, it is assumed that encounters in the blood between T_e cells and C cells are modeled by the law of mass action, i.e., the effect of these cells on each other is based on the sizes of the two cell populations, and there is no saturating effect in the blood circulation system. The lower case parameters (S_n, an, etc.) in the above equations are all constants, as is C_max [2].

Figure 1: Cell population diagram. In this figure, solid curves with arrows represent source terms (such as proliferation and activation), dotted curves without arrows represent interaction terms (such as contact between cells), and dashed curves with arrows represent loss of cells (such as death). Curves with arrows signify movement into or out of a population, while curves without arrows only signify interactions between populations [14].

Here, we use the same values for the parameters as in [14] and appear in table 1.

Table 1: Parameter Values.

As with ODE, the dynamic status of system (1-3) can be studied. First, we should find the fixed points of system (1-3) and then their stability should be analyzed. In system (1-3) we reach the problem parameters to 8 using standard rescaling (non-dimensionalization) [ 14]. T_n is rescaled with factor d_n/S_n, T_e is rescaled with factor γ_c/d_n, and C is rescaled with factor γ_e/d_n and t is rescaled with factor d_n. With these rescaling system (1-3) yields

Where the new factors are defined as

For finding the fixed points considering to C=0 (no cancer cell), Eq. (4) implies that T_n=1 and Eq. (5) implies that T_e=0. There are no any other fixed points for C=0. So, we have just P₁= (1, 0, 0). To find the other fixed points, for the case C≠0, from Eqs. 4-6 we have

The third term on the right-hand side of Eq. (7) is logarithmic in C and hence, increases as C increases. For the factors of the fourth term to be zero, C must be negative (for the clinically feasible ranges of parameters in table 1). Thus, this fourth term is a rational function which decreases as C increases (when C is positive), and so there is at most one value of C on which makes the right-hand side of Eq. (7) to be zero. Therefore, the second fixed point will be

and for this we need

Note that again, for biological purpose the populations of T_n, T_e and C should not be negative. Therefore, we have at most one real fixed point, rather than P₁ for this system.

To determine the stability analysis of the cell populations near the fixed points, according to the Matignon theorem [18], the fixed points of FDE system (4-5) is asymptotically stable if and only if the eigen values of related linearized system satisfied in

. However, since the minimum value of α ∈ (0,1) that we chose here is closed to one, namely α = 95.0, the stability domain for the eigen values of the linearized system will be (almost) left hand side of R² coordinates, i.e. the place that the real part of eigen values are negative. This means that we may use the same theorems as in ODE systems for the stability analysis of our FDE system.

Therefore, for the stability analysis of the fixed points, we need to determine the linearization of the system (4-6). This linearization yields

By substituting P ₁= (1, 0, 0) in this matrix we get

By easy calculation, the eigen values of this Jacobian matrix will be λ₁= −1, λ₂= −ξ₅ and λ₃= −ξ₆ − ξ₈. It is clear that λ₁, λ₂ and λ₃ all have negative (real) sign Therefore, P₁ is stable. For the second fixed point

(if exists), we may use a similar analysis as above. However, due to the long and complicated calculations, such as the one has done by Routh test [19], presents a kind of non-algebraic system that cannot be solved analytically [13].

As an alternative, we can calculate the eigen values for a wide variety of possible parameter values, as has done in [14], by systematically sampling through the ranges listed in table 1. Using this method to find the eigen values of the matrix (8), we may use different sampling intervals for

and

. For example, using the interval (0, 5000) for populations ,

and (1, 400000) for

all of the eigen values were either negative real or were imaginary with all of the real parts bounded above by -1.000076. The large ranges of

and

used to calculate the eigen values give reasonable confidence that the eigen values have negative real parts for this fixed point. Hence, we assume the second fixed point to be asymptotically stable.

Discretization and Numerical Solution in Fractional Mode

As we discussed above, linear stability analysis of system (1-3) or (4-6), around its fixed points, were similar to that of its ODE counterpart. However, to solve FDE system (1-3) first we need to discretize it. Among the several discretization methods that are available for the fractional derivative

, we used the one that have generated by Grunwald- Letnikov [16,17]. In this method

is approximated by

where, l is the step size and [t] is the integer part of t. Using this method for system (1-3),

is replaced by

where t_n = nl and

is Grunwald-Letnikov coefficients defined by

We may calculate

with the following recursive formula too.

Now, Using Eq.9, system (1-3) discretize as follows.

By simple calculations, system (11-13) yields the following recursive formulas.

The numerical results carried out by using MATLAB™ software for two sets of parameters given in table 1 (patients A and B) with initial values (T_n)₀=1510, (T_e)₀=10 and C₀=10000. These results illustrated in figures 2- 9 for both patients, whom we label A and B for different values of derivative order α. Each set of parameters defines an individual (hypothetical) patient for whom we determine an optimal dosing strategy. As we can see in figures 2- 5, the graphs shown by dashpoints (black color) present tumor concentration, C, for patients A and B with derivative orders α=1 and α=0.95. If we look at these graphs more precisely, we claim that the ones belong to the fractional order derivatives (α=0.95) behaves more naturally than the one belong to classical order (α=1). For example, comparing the C carves in figures 4 and 5, the one belong to classical order goes from initial value 10000 to the minimum point 8500 (Figure 4). While, the one belong to the fractional order goes to the minimum point 8000 and then growth up by a natural slop. It seems that the model with fractional order is closer to the nature of the body defense system. However, this claim should be verified by more precise clinical data.

Figure 2: Numerical results for cancer cell population (C) (patient A) found from ODE system for classical derivative order α=1.

Figure 3: Numerical results for cancer cell population (C) (patient A) found from FDE system for fractional derivative order α=0.95.

Figure 4: Numerical results for cancer cell population (C) (patient B) found from ODE system for classical derivative order α=1.

Figure 5: Numerical results for cancer cell population (C) (patient B) found from FDE system for fractional derivative order α=0.95.

Figure 6: Numerical results for naive T cell population (T_n) (patient B) found from ODE system for classical derivative order α=1.

Figure 7: Numerical results for naive T cell population (T_n) (patient B) found from FDE system for fractional derivative order α=0.95.

Figure 8: Numerical results for effector T cell population (T_e) (patient B) found from ODE system for classical derivative order α=1.

Figure 9: Numerical results for effector T cell population (T_e) (patient B) found from FDE system for fractional derivative order α=0.95.

FDE Model with Drug Treatment

In this section, the ODE model presented by Moore and Li for brain blood cancer [14] is derived using FDE. So, if we consider the same system (1-3) with three cells populations along with a chemotherapy treatment describing the growth, death, and interactions of each cells, then such this system can be formulated by the means of FDE as follows.

In this system (0), (0) and C(0) are known initial values and time dependent drug efficacies are incorporated using u₁(t) and u₂(t). Setting u ₁(t) ≡ 0 and u₂(t) ≡ 1 in the Eqs. (17-19) would give the same model described for the dynamics of the disease without treatment. All of the parameter values in these equations are assumed to be positive. Again, the structure of the equations guarantees non negative solutions for the state variables, (t), (t) and C(t). The negative terms in the above equations represent losses from the cell populations while the positive terms are source terms for the cell populations [2].

The effect of the targeted drug represents by the control function u₁(t), which slows the production of cancer cells. We assume this drug affects only cancer cells and not the other cells, so u₁(t) appears only in Eq. (19). The u₂(t) term uses to incorporate treatment by a broad chemotherapy, such as cytarabine or hydroxyurea or a combination of such drugs, which is cytotoxic to all three-cell populations. Thus, u ₂ appears in all three state equations as a coefficient in the natural attrition terms. Values of u₂>1 correspond to treatment with a cytotoxic drug [2]. In this case, cells T_n, T_e and C decreases with constant factors d _c, d_e and d_n, respectively. But when the patient is treated by drug, considering−u ₂(t)d_nT_n, −u₂(t)d_eT_e and −u₂(t)d_c then T_n, T_e and C are decreased more. CML cells are more decrease in patient B than A considering d_c is higher than d_e and d_n. In this case, we are minimizing the value of u₁ and u₂ subject to the equations of system (17-19). That is ((u₁ ,u₂ ) ∈U),

Where U={(u₁(t),u₂(t)) s.t. m_i<u_i(t)<M_i,u_i lebesque measurable i=1,2, t ∈[0,t_f]}.

In the objective functional, we minimize the total cancer cell population over the interval [0, t_f] through the first term in the integrand, and at the final time through a salvage term B₃(t_f). We also minimize the systemic costs to the body of the two drugs u₁ and u₂ . As in [16,17], it is expected that the effects of the drugs are non-linear, and we choose quadratic cost terms

to reflect these effects. The coefficients B₁ and B₂ are weight constants on the controls, and include a measure of toxicity of the drugs to the body. We note that the higher the weight the greater will be the toxicity. The salvage term B₃(t_f) term was not present, the controls could taper off earlier, and allow a rise in cancer cell count at the end of the treatment period. The salvage term −B₄(t_f) included to penalize for low values of T_n, since this affects the patient’s ability to fight off other diseases. The coefficients B₃ and B₄ allow the salvage terms to be weighted differently from each other and the integral terms. (The coefficients B₁, B₂, B₃ and B₄ are all positive.) The lower bounds for u₁ and u₂ correspond to no therapy. For u₁ this lower bound is m₁=0, and for u₂ the lower bound is m₂=1. We suppose M ₁<1, as M₁=1 would correspond to no new cancer cells. The upper bound M₂ is greater than 1 and is determined by the parameters d_c, d_e and d_n [2] in such a way that is obtained by

The following hypothesis is used for optimizing u₁(t) and u₂(t).

Theorem 1 (Characterization of the Optimal Control)

Suppose an optimal control

and the solutions of system (17-19) that minimize the function

are given. Then there exist adjoint variables λ_i for i=1,2,3 and 4 that satisfy to the following system

Here, λ₁(t_f) = −B₄, λ₂(t_f) = −B₄ and λ₃(t_f) = −B₄. Moreover, the optimal

is given by

Since, the state and adjoin solutions are a priori L∞-bounded, the right-hand side of the state and adjoint equations become Lipschitz in those solutions [14]. This Lipschitz property guarantees that the solution of the optimality system is unique if the final time is sufficiently small. We can see the paper by Fister et al. [20] for a uniqueness proof using Lipschitz properties. The uniqueness of the solutions of the optimality system implies the uniqueness of the optimal control pair.

Now for solving optimization problem (20), we should first solve system (21-23) with some initial values of T_n, T_e and C. Here, to be consistence with other results in (16,17), we start with T_n=1510, T_e=10 and C=10000. Then, by finding the value

and plugging into the system (17-19), instead of u₁(t) and u₂(t), we are ready to solve 12 this system with the same starting point T_n, T_e and C as above. Similar discretization method that we have done for FDE system (1-3) can be applied here for system (17-19) to get

Now by simple calculation on (24-26), we get the following recursive formula.

By solving this system for some customary time, say t∈[0,1], we arrive at the new point T_n, T_e and C. Then, this new set of values will serve as new starting point with initial values λ₁(t_f)=−B₄, λ₂(t_f)=0 and λ ₃(t_f)=−B₃ for solving system (21-23), in the next iteration, to find a new optimal value of

. These iterations will continue up to the time t=250 (days). Indeed, using MATLAB™ to solve these two joint systems (27-29) and (21-23), as the solution of optimal problem (20), the results are illustrated in figures 2- 9 for different values of fractional derivative 0.90<α≤1. The results of patients A and B for α=1 (classical ODE) are shown in figures 2, 4, 6 and 8. These results are the same as in reference [2]. Figure 2 compares the behavior of CML cell population for one control (u₁) as well as two controls (u₁ and u₂). For a more aggressive case of cancer (as for patient A) the CML cell count is significantly lower after treatment with combination therapy as compared to treatment by targeted therapy alone. For a less aggressive case (patient B) in figure 4, it may suffice to treat with one drug alone as the outcome is similar for single drug therapy and combination therapy. However, the CML plot of patient A is shown in figure 3 for α=0.95, which is different from figure 2.

As we can see in figure 2, with just one treatment u₁, the CML cells increase to their lower bound 0.8 after 250 days. While in figure 3, which shows the results for our FDE model, these cells increase to the value 0.4. With the same comparisons in figure 2, when we have both treatments u₁ and u₂ the increasing slop of CML is lower than the slop for this cells in figure 3. We claim that the results found in our FDE model is closer to the nature of the drug treatment than the results found by ODE contra parts. However, this claim should be verified by the clinical treatment data.

The CML plot of patient B is shown in figure 5 for α=0.95, which is again different from the plot in figure 4 for the same patient B using the classical ODE. In figure 6, there is not much difference in T_n evolution over time with or without drug therapy. Under targeted therapy, immune response is not compromised due to drug dosing. The T_n plot of patient B for fractional derivative α=0.95 is shown in figure 7, which is again different from the classical ODE plotted results in figure 6. These different are clear from the slops of the curves T_n in both figures.

Finally, we see in figure 8 for patient B, whose cancer is of a less aggressive nature, that T_e cell response is enhanced when treated only with u₁. Using our FDE model, the similar results for patient B plotted in figure 9 with fractional derivative α=0.95. As we can see in figures 8 and 9, there is no significant different between the results whenever we have no treatments or both treatments u₁ and u₂ are on. However, for one drug treatment u₁, in the case of classical ODE model the maximum of T_e will be 24, while in FDE model this maximum will be 38. After these maximum, in both cases, T_e will converge to zero.

We claim that the results of figures 3, 5, 7 and 9 are more consistent with the nature of drug therapy. We emphasis, one reason for the accuracy of our FDE model is the non-local property of fractional derivative. This means that the next state of a system not only depends upon its current state but also upon its historical states starting from the initial time. To see this, pay attention to the summation terms in right hand side of system (27-29). However, as we have said above, this claim should be verified by more clinical treatment data.

Conclusion

In this article, we have studied a mathematical model with fractional order derivatives as a dynamic system for presenting the transaction between body immunology and drug variable. We have introduced a three cells population model describing the interaction between the CML and the naive T cells together with the effector T cell population, without any treatment, in the form of FDE. As we have seen, the local stability analyses of the fixed points for this FDE system were the same as its counterpart ODE system. These analyses were agreed with numerical results of discretized FDE system using Grunwald-Letnikov method. As we have expected the tumor cell population were increasing up to its maximum values by a positive initial value. Hence, a more reliable FDE system with chemotherapy treatment was considered. In order to find the best amount of medicine on which the tumor cell population, the naive T cells and the effector T cell population, decreasing we conducted a optimal control. The drug optimized dose is resulted from targeted therapy and broad cytotoxic therapy. We could adapt the same existence and characteristic optimal control theorems as in ODE systems for our FDE system. We claim, due to the non-local property of FDE, the results found by this system were more accurate as we compare to the results found by counterpart ODE models. Of course, this claim should be verified by more clinical treatment data.

Here, we should emphasis that by choosing the small values of α∈[0,1] , we will encounter to the larger amount of error in calculations. In FDE models that we have introduced here, experimentally, we have found that the best value of α for the best results is 0.95.

Acknowledgement

This work was supported by Health Policy Research Center, Shiraz University of Medical Sciences

References

Faderl S, Talpaz M, Estrov Z, O’Brien S, Kurzrock R, et al. (1999)

Nanda S, Moore H, Lenhart S (2007)

Deininger MW, Goldman JM, Melo JV (2000)

Deininger MW, O’Brien SG, Ford JM, Druker BJ (2003)

Deininger MW, Druker BJ (2003)

Van Rhee F, Szydlo RM, Hermans J, Devergie A, Frassoni F, et al. (1997)

Baccarani M, Martinelli G, Rosti G, Trabacchi E, Testoni N, et al. (2004)

Kano Y, Akutsu M, Tsunoda S, Mano H, Sato Y, et al. (2001)

Topaly JW, Zeller J, Fruehauf S (2001) Synergistic activity of the new ANL-specific tyrosine kinase inhibitor STI571 and chemotherapeutic drugs on BCR-ABL-positive chronic myelogenous leukemia cells. Leukemia 15: 342.

Thiesing JT, Ohno-Jones S, Kolibaba KS, Druker BJ (2000)

Fokas AS, Keller JB, Clarkson BD (1991)

Adimy M, Crauste F, Ruan S (2005)

Afenya EK, Bentil DE (1998)

Moore H, Li NK (2004)

Einstein A (1998) Mount Sinai School of Medicine. Annual Publication of the American Institute of Physics.

Podlubny I (1999)

Oldham KB, Spanier J (1974) The fractional calculus, New York: Academic Press Inc, 234.

Matignon D (1996) Stability results of fractional differential equations with applications to control processing, in: IMACS, IEEE-SMC, Lille, France, 963.

Borrelli R, Coleman C (1998) Differential Equations: A Modeling Perspective. Wiley, New York.

Fister KR, Lenhart S, McNally JS (1998)