Search tips
Search criteria 


Logo of springeropenLink to Publisher's site
Journal of Pharmacokinetics and Pharmacodynamics
J Pharmacokinet Pharmacodyn. 2010 June; 37(3): 257–276.
Published online 2010 May 9. doi:  10.1007/s10928-010-9159-z
PMCID: PMC2889283

Fractional dynamics pharmacokinetics–pharmacodynamic models


While an increasing number of fractional order integrals and differential equations applications have been reported in the physics, signal processing, engineering and bioengineering literatures, little attention has been paid to this class of models in the pharmacokinetics–pharmacodynamic (PKPD) literature. One of the reasons is computational: while the analytical solution of fractional differential equations is available in special cases, it this turns out that even the simplest PKPD models that can be constructed using fractional calculus do not allow an analytical solution. In this paper, we first introduce new families of PKPD models incorporating fractional order integrals and differential equations, and, second, exemplify and investigate their qualitative behavior. The families represent extensions of frequently used PK link and PD direct and indirect action models, using the tools of fractional calculus. In addition the PD models can be a function of a variable, the active drug, which can smoothly transition from concentration to exposure, to hyper-exposure, according to a fractional integral transformation. To investigate the behavior of the models we propose, we implement numerical algorithms for fractional integration and for the numerical solution of a system of fractional differential equations. For simplicity, in our investigation we concentrate on the pharmacodynamic side of the models, assuming standard (integer order) pharmacokinetics.

Keywords: Link model, Direct action, Indirect action, Fractional calculus, Numerical methods


Since Leibniz’s prophecy in the late seventeenth century to a few decades ago, expressions involving fractional derivatives, integrals and differential equations were mostly restricted to the realm of mathematics. The first modern examples of applications can be found in the classic papers by Caputo [1] and Caputo and Mainardi [2] dealing with the modeling of viscoelastic materials, but it is only in recent years that it has turned out that many phenomena can be described successfully by models using fractional calculus. In physics fractional derivatives and integrals have been applied to fractional modifications of the commonly used diffusion and Fokker–Planck equations, to describe sub-diffusive (slower relaxation) processes as well as super-diffusion [3]. Other examples are of applications are in diffusion processes [4], signal processing [5], diffusion problems [6]. More recent applications are mainly in physics: finite element implementation of viscoelastic models [7], mechanical systems subject to damping [8], relaxation and reaction kinetics of polymers [9], relaxation in filled polymer networks [10], viscoelastic materials [11], although there are recent applications in splines and wavelets [12, 13], control theory [14, 15], and biology [16] (bacterial chemotaxis). Surveys with collections of applications can also be found in Nonnenmacher and Metzler [17], and Podlubny [18]. The suggestion of possible PK applications is reported in [19, 20], and generalized in [21], while, to the best of our knowledge, there are no applications to pharmacodynamics.

The reasons for the lack of application to pharmacodynamics are mainly two: first and foremost, PKPD models incorporating fractional calculus have not been proposed; second, it turns out that one needs to provide special algorithms to deal with such models. In this paper we propose new families of PKPD models incorporating fractional calculus, and investigate their behavior using purposely written algorithms. The main purpose of this paper is to attract the attention of the field into these possibly interesting families of PKPD models.

The structure of this paper is as follows. In the next section we briefly recall the mathematical foundations of fractional calculus, and summarize the algorithms for the numerical methods used in the solution of problems described by fractional-order integrals, and differential equations. We then provide illustrations of the application of fractional order integrals, and systems of fractional order differential equations PKPD models. In particular, we will focus on link models, direct concentration–response relationships describing models relating drug concentration (in plasma or biophase) to a pharmacodynamic effect, and indirect mechanism of action models that involve primarily drug induced modulation of endogenous factors. A final discussion comments on the results, the general families of models we introduce, and the challenges associated with the application of the models to real data. A first appendix provides (simplified) code for the implementation of the numerical methods used in the paper, while a second appendix describes some numerical properties of fractional integration which might help in the interpretation of applications to PD data.

Fractional integrals and derivatives

The calculus of fractional integrals and derivatives is almost as old as calculus itself going back as early as 1695, to a correspondence between Gottfried von Leibnitz and Guillaume de l’Hôpital. About 300 years had to pass before what is now known as fractional calculus was slowly accepted as a practical instrument. A brief history of the development of fractional calculus can be found in Miller and Ross [22], while a survey of many emerging applications of the fractional calculus in areas of science and engineering can be found in the text by Podlubny [23], here we give a brief overview.

Although mathematical modelers dealing with dynamical systems are very familiar with derivatives of integer order, equation M1 and their inverse operation, integrations, they are generally much less so with fractional-order derivatives, for example equation M2 One way to formally introduce fractional derivatives proceeds from the repeated differentiation of an integer power p:

equation M3

For p a real number, repeated differentiation gives

equation M4

where the gamma functions [24] replace the factorials. The gamma function allows for a generalization to a non-integer order of differentiation, with α real, equation M5

equation M6

The extension defined by Eq. 3 corresponds to the Riemann–Liouville derivative [4, 22].

A more elegant and general way to introduce fractional derivatives uses the fact that the m-th derivative is an operation inverse to m-fold repeated integration. Basic to the definition is the integral identity

equation M7

Clearly, the equality is satisfied at x = a, and it is not difficult to see iteratively that the derivatives of both sides of the equality are equal. A generalization of the expression allows the definition of a fractional integral (FI) of arbitrary order via

equation M8

where again the gamma function is replacing the factorial. In this paper we are concerned with fractional time derivatives, and we take the lower limit in Eq. 5 to be zero. For this reason in the following we will drop the subscript a in the definition of the operators we consider, and use t, instead of x, to indicate the independent variable time. Starting from Eq. 5, one can construct several definitions for fractional differentiation. The fractional differential operator equation M9 is defined by

equation M10

where m is the smallest integer greater than α, equation M11 (equation M12 integer) is the classical differential operator, and f(t) is required to be continuous and α-times differentiable in t. The operator equation M13is named after Caputo [1], who was among the first to use it in applications and to study some of its properties. It can be shown that the Caputo differential operator is a linear operator, i.e. that for arbitrary constants a and b:

equation M14

which is important in PK and PD applications to, e.g., preserve linearity in respect to inputs; that it commutes:

equation M15

which is important to compute partial derivatives in respect to parameters for model fitting, or for future applications involving mixed effect models [25]; and that, as the standard differential operator, it possesses the desirable property that:

equation M16

for any constant c, e.g. the Caputo derivative for steady-state values is zero, as it should for physical realism.

Having defined equation M17 we can now turn to fractional differential equations (FDE), and systems of FDE, which will be used to specify PKPD models and will need to be solved over an interval [0,t], in accordance with appropriate initial conditions. A FDE of the Caputo type has the form

equation M18

where y(t) is a vector of dependent state variables, and f(t,y(t)) a, dimensionally conforming, vector valued function, satisfying a set of (possibly inhomogeneous) initial conditions for the y(t) and its derivatives

equation M19

equation M20 where equation M21 is the ceiling function, giving the smallest integer greater than or equal to its argument. It turns out that under some very weak conditions placed on the function f of the right-hand side of Eq. 10, a unique solution to Eqs. 10 and 11 does exist [26].

A typical feature of differential equations (both classical and fractional) is the need to specify additional conditions in order to produce a unique solution. For the case of Caputo FDE, these additional conditions are just the initial conditions listed in Eq. 11 which are similarly required by classical ODEs. In contrast, for Riemann–Liouville FDE, these additional conditions constitute certain fractional derivatives (and/or integrals) of the unknown solution at the initial point t = 0 [27], which are functions of t. These initial conditions are not physical; furthermore, it is not clear how such quantities are to be measured from experiment, say, so that they can be appropriately assigned in an analysis [22]. If for no other reason, the need to solve FDE is justification enough for choosing Caputo’s definition for fractional differentiation over the more commonly used (at least in mathematical analysis) definition of Liouville and Riemann, and this is the operator that we choose to use in the following.

Numerical approximation of caputo-type derivatives

One of the main factors limiting the applications of fractional order integrals, and differential equations to the modeling of real data is the fact that analytic solutions to equations of the form (5) and (10)–(11) are know only in simple cases, or for specific values of α, see for example [28, 29]. As a consequence one is forced to use numerical approximations to obtain the desired solutions. However, to date, algorithms to do so have been developed only to a rather limited extent, and none is available either commercially or as free software. To investigate the applications of FI and FDE to PK/PD models we implemented numerical methods reported as algorithms by [30, 31]. As indicated above, the natural concept for the fractional integral to be used in connection with Caputo differential equations is the Riemann–Liouville FI described in Eq. 5, and, as we will see, Riemann–Liouville integrals will have applications in PKPD models. We therefore implemented an algorithm for the solution of FI problems, which employs a product integration technique based on the trapezoidal quadrature rule. That is, one replaces the given function f(t) by a piecewise linear interpolant, and then we calculate the resulting integral exactly. The algorithm is also part of the numerical solution of FDE, which is of the Predict–Evaluate–Correct–Evaluate type and is reported in simplified form the Appendix 1 as R code [32]. Appendix 2 shows some properties of the integration algorithm and fractional integrals, and we refer to [30, 31] for further details. Other numerical algorithms exist that solve FDE, e.g., [23, 33], but they focus on solving Riemann–Liouville FDE and usually restrict the class of FDE to be linear with homogeneous initial conditions, the algorithm reported in the appendix solves non-linear Caputo differential equations with inhomogeneous initial conditions, if required. The computations reported in this article are done using a tolerance of 0.001 (see Appendix 1).

Fractional Integrals, equation M22 and Caputo differential equations, equation M23 are the operators which are used to generalize the standard PKPD models, as described below.

Pharmacokinetics–pharmacodynamics models

A pharmacodynamic model is traditionally divided in two parts: a link model and a transduction model (see [34] for a detailed description and review of pharmacodynamic models). The link model describes the distribution of drug from an observed compartment (e.g. plasma) into a biophase (the “effect compartment”). Drug in the biophase induces a pharmacodynamic response (the “effect”), which is described by a transduction model. We will consider two types of transduction models: direct-response and indirect-response. In the case of direct response models, the effect is represented by a non-linear memory-less function of drug concentration in the effect site (Ce, the model was introduced by [35]). For indirect action models the action of drug is also represented by a non-linear memory-less function of Ce, in this case acting on the rate of production or elimination of an endogenous substance (the effect) [36]. See [34] for a detailed description and review of pharmacodynamic models.

In the following we assume that drug kinetics in plasma are represented by standard (integer order) differential equation models, or their analytical solutions.

Link model

The model describing the transfer of drug to the effect compartment is, it its mechanistic interpretation, part of the pharmacokinetics of the drug. However, in practise, and because drug concentration in the effect compartment is not observed, it is often unclear if the link model is actually associated with a site of action, or if it represent some pharmacodynamic related delay, or a combination of the two. We start by representing drug concentration in the effect compartment by the (Caputo) fractional differential equation:

equation M24

where equation M25represents the drug input into the central compartment and equation M26 is the elimination rate constant from the effect compartment. In general the solution of Eq. 12 takes the form of the convolution:

equation M27

where equation M28 and equation M29is the Mittag–Leffler function [37]:

equation M30

The Mittag–Leffler function is related to the monoexponential by the relationship equation M31 and it represents the unit-response function for fractional kinetics (see [21] for general expressions of the unit-response functions for arbitrary fractional compartmental structures).

Note that when Ce0 = 0,equation M32 as for the standard effect compartment model, but it also holds equation M33 a fact that might induce a-posteriori identifiability problems.

In the following example equation M34 represents drug concentration in the central compartment following, e.g., an oral administration of a dose Dose; in the example Dose = 1, the initial condition in the effect compartment is Ce0 = 0, and parameters values are V = 1, ka = 5.0, kel = 0.5, and keo = 1. The upper panel of Fig. 1 shows concentration in the effect compartment for α = 1 (thick solid line, computed using the analytic solution to model (12)), and the solution obtained using the predictor–corrector algorithm (superimposed thin solid line, indistinguishable from the analytic solution). The additional dashed lines show the solution to model (12) with α = 0.9, 0.75, 0.5, and 0.25, computed using the predictor–corrector algorithm. One can notice the delay in the peak concentration of Ce(t) (time to peak concentration, tpeak, is at 1.63 (unit time, ut)) versus Cp (tpeak = 0.50 ut), delay that shortens as α decreases (tpeak = 1.47, 1.28, 0.97, 0.69 ut). In addition the decline is slower than exponential. The lower panel of Fig. 1 shows Ce(t) profiles for α = 1, 0.25, 0.01, 0.001, 0.0001, and 0.00001. Following the profiles for decreasing values of α, note the change of Ce(t) from initial convexity to concavity, and its progressive approach to the profile of drug concentration in plasma (Cp(t), thick dashed line), which is exactly reached for α = 0.0, when equation M35 the Dirac delta-function. A feature of the family of curves corresponding to α < 0.01 is that tpeak in the effect comportment is earlier than in the central compartment: tpeak is approximately equal to 0.48, 0.32, 0.45, for α = 0.01, 0.001, 0.0001, respectively, tpeak = 0.50 for α = 0.00001.

Fig. 1
Numerical solutions to the link model equation (12) versus time (t). equation M36, Ce (0) = Ce0, f(t)=equation M37, for different values of α (see text for parameters values). Upper panel (thick solid line): analytic solution for α = 1, ...

Direct action models

In the standard direct action model the effect at time t, equation M38 is expressed by an arbitrary (memory-less) function of drug concentration in the effect site at time equation M39 however to generate a wider class of relationships, we assume that the effect at time t is related to the Riemann–Liouville fractional integral of Ce(t), a quantity that we call acting drug:

equation M40

For a direct action model the effect is simply related to A(t) by a non-linear memory-less function G:

equation M41

Equation 16 for β = 0 represents a standard direct action model, the acting drug is drug concentration in the effect site: equation M42 for β = 1 it represents an effect related to the exposure,equation M43, where equation M44; for equation M45 the effect is related to a fractional integral that is in between Ce(t) and exposure, i.e., it represents effects depending on hypo-exposure to drug, for equation M46 we have the cases of hyper-exposure dependencies.

The upper panel of Fig. 2 shows equation M47 where Ce(t) is computed as before (α = 0, i.e. the effect compartment has first order kinetics) for different values of β, using the numerical approximation to Riemann–Liouville integrals. Note that equation M48 for equation M49, while for β = 1 equation M50.

Fig. 2
Numerical solutions to the direct action model equation (16) versus time (t). equation M51, Ce(t) computed as in Fig. 1. Upper panel: numerical solutions to equation M52 for β = 0 (equation M53, i.e. drug concentration in the effect compartment), β = 0.1, ...

The lower panel of Fig. 2 shows the corresponding effect profiles when:

equation M56

and E0 = 0.1, Emax = 1, Ce50 = 0.15. Note how for β = 0, 0.1, 0.3, 0.6, and 0.9 the effect return to baseline, equation M57, but the return is progressively delayed. For β = 1, the effect, reaches the maximum equation M58, and does not return to baseline.

The upper panel of Fig. 3 shows equation M59 for β = 1, 1.25, 1.5, 1.75, 2.0, i.e., cases in which the effect is either proportional or hyper-proportional to exposure to drug. In the upper panel, note how the Riemann–Liouville fractional integration introduces a delay in the increase of the active drug for β > 1, and that while for β = 1, A(t) asymptotes to Dose, that isequation M60, for β > 1 we haveequation M61. The lower panel of Fig. 3, shows the effect profiles computed according to (17), where now, for β > 1, equation M62.

Fig. 3
Numerical solutions to the direct action model equation (16) versus time (t). Upper panel: numerical solutions to equation M63, for β = 1, exposure to drug concentration in the effect compartment, for β = 1.25, 1.5, ...

Indirect action models

Dyneka et al. proposed four indirect response models [36] that are based on drug (Cp(t) or Ce(t)) effects that either stimulate or inhibit the production or loss of an effect (or its mediator). The models take the form:

equation M66

where kin and kout are the production and elimination rate of Y(t), and, G(Ce(t)) positive drug stimulates or inhibits the production rate depending on the sign in the equation; or

equation M67

where drug stimulates or inhibits the elimination rate of Y(t). Fractional dynamics can be immediately introduced by substituting equation M68, or equation M69, for Ce(t) in Eq. 18. However, an additional alternative class of models is generated if we make the dynamics of Y(t) to be represented by a fractional differential equation. We first consider an inhibition of production rate:

equation M70

when kin = 1, kout = 1, equation M71, Ce(t) is as in the previous examples. The upper panel of Fig. 4 shows the effects profiles corresponding to model (20), for equation M721. The main features are the quicker onset and slower decay when equation M73, dashed lines, in respect to the standard model equation M74, solid line. The lower panel of Fig. 4 shows effects profiles for equation M75. Now when equation M76the onset of the effect is delayed, and damped oscillations are introduced on the way to the return to baseline kin/kout = 1.

Fig. 4
Numerical solutions to the indirect action model equation (20) versus time (t). equation M77, equation M78, equation M79, Ce(t) computed as in Fig. 1 (see text). Upper panel (solid line): Y(t) corresponding to the standard indirect action model, δ = 1; dashed ...

Using a similar layout, Fig. 5 shows the effect profiles corresponding to an inhibition of elimination rate:

equation M80

where all the parameters are as in Fig. 4, but equation M81 with γ = 2. Note how in Fig. 4 and and55 the Y(t) can take negative values, indicating that, depending on the application, γ will need to be constrained to satisfy the physical properties of the response.

Fig. 5
Numerical solutions to the indirect action model equation (21) versus time (t). equation M82, equation M83, equation M84, Ce(t) computed as in Fig. 1 (see text). Upper panel (solid line): Y(t) for δ = 1; dashed lines: Y(t) for 0 < δ < 1. ...

The main observations that can be derived from the examples so far reported as well as additional extensions are discussed below.


Fractional calculus has only recently came to the attention to pharmacokinetics modelers. The use of the Mittag–Leffner function, instead of a monoexponential, to describe single-compartment kinetics following bolus input was suggested in [19], [20] extends the approach to commensurate1 two compartmental models, while [21] points out that certain sums of single- and two-parameters Mittag–Leffler functions turn out to be the response functions corresponding to rather arbitrary, commensurate and non-commensurate, compartmental systems [21], those sums, and their convolution with arbitrary inputs, can be used to represent a very large class of fractional pharmacokinetics models.

This paper turns the attention to pharmacodynamics, and in so doing it introduces a new approach to pharmacodynamics modeling using fractional integrals and fractional differential equations. The approach not only extends well-known direct- and indirect-action models, but also generates novel families of PD models, which are reported, for easy reference and summary, in Table 1. The models we propose offer mechanistic insight very similar to the ordinary ones, and in particular models parameters, other than the fractional orders, allow the same interpretation.

Table 1
Fractional calculus based link, active drug, and pharmacodynamics models

The examples reported above are specific instances of these families: Fig. 1 shows the link models when α varies; Fig. 2 and and33 direct action models for α = 1, β variable, and δ = 1; Fig. 4 and and55 indirect action models for α = 1, β = 1, and δ variable. We now comment further on the models we introduce.

Link models

In addition to possible additional variety of kinetics for the effect compartment introduced by the use of Eq. 12, the observed behaviour of the link models (Fig. 1) suggest a possible application to situations in which the observed peak effect is observed earlier than the peak concentration. The example of nicotine effect on heart rate represents one of these possible applications [38].

Direct action models

The main novelty demonstrated by Fig. 2, is the introduction of what we call “active drug”, the fractional integral of the concentration in the effect compartment. The active drug serves as independent variable for the pharmacodynamic model, and this, as shown by Eq. 16 and Fig. 2, allows the definition of a large family of models. The models smoothly range from a standard dependency on Ce(t) (β = 0), to dependency on hypo-exposure (0 < β < 1), to a standard exposure model (β = 1), and dependency on hyper-exposure (β > 1). The family generated by varying β potentially greatly expands the current possibilities, limited to β = 0 or 1, in addition this is achieved using the simplest tool of fractional calculus (fractional integration). Hypo-exposure models are particularly interesting because they establish a relationship with integrals of drug concentrations, but they are reversible: for 0 < β < 1 eventually the effect returns to baseline when drug is eliminated from the body, thus this feature removes one of the main limitations of models depending on exposure. Hyper-exposure models are, as well as exposure models, not reversible, but they might have applications, for example in toxicity studies. (Appendix 2 further comments on the properties of fractional integration in its relationship to the concept of “active drug” and hypo-/hyper-exposure models.)

Indirect action models

The use of fractional dynamics as in models (20) and (21), allows to smoothly transition from models with exponential or sub-exponential tails (equation M96), to models showing damped-oscillations (1 < α < 2). The latter could offer interesting applications to effects exhibiting rebound and (damped) oscillatory return to baseline, possibly avoiding the formulation or more complicated models describing such behaviours.

Setting α = 1, but using active drug, equation M97, instead of Ce(t), introduces a family of models similar to the direct action model described by Eq. 16, one that was not considered in the examples reported above.

General considerations

As one reviewer pointed out “a great benefit for the methods proposed is increase flexibility of response-time profiles with only a slight expenditure in the number of parameters. That is, the FDE approach appears to lend itself to parsimony. The value of this flexibility is appealing, especially if one considers PKPD models typically postulated as simplifications of complex biological process.” The reviewer also pointed out that it would be useful to relate the FDE approach to a mechanistic model. For example, the use of a Mittag–Leffner function, instead of a mono-exponential, can be conceptually related to more complicated diffusion processes than traditionally assumed [3]. For the PD models we introduce, the same argument can be made for the FDE link model, while for the FDE extension of direct/indirect-action models one might associate the FDE representations to more complicated binding or transduction processes. It is also possible that the models transitioning into damped-oscillations (1 < α < 2) could be used to as lamped approximation to feed-back type processes. Fractional exposure models have possibly the simplest explanation, representing a very natural generalization of traditional exposure models. Mechanistically a fractional exposure model might be generated, for example, in presence of modulation of excretory transport [39].

Open problems and questions

The main problem to solve to apply the models we present in this paper to real data is the development of flexible algorithms that can be used in conjunction with an estimation program such as NONMEM [40] or NLME [41] for model identification. In particular the algorithm used here [30] will need to be modified to allow for non-equispaced observations as well as to allow for multiple fractional components. We are collaborating with Prof. Kai Diethelme and Neville Ford to develop more robust and flexible FI and FDE numerical solvers.

As the examples reported in this paper demonstrate, the behaviour of FI and especially FDE are quite sensitive to small variations in α, and this sensitivity might pose problems when estimating α as well as the other parameters of a PKPD model. The usual identifiability issues [42], and the use of constraints will need to be investigated as well.

It is quite obvious, in reference to Table 1, that one could devise more complicated examples than the ones shown and discussed above: PKPD models are well defined for any combination of α, β, and δ. (For example, for a direct action model setting α,β = 1, but δ variable, introduces yet another family in which the effect is related to a fractional integral of the function equation M98) While the flexibility provided by the possible choices in α, β, and δ should offer a rich and interesting family of models, practical problems could arise simply due to the large number of models that would need to be considered. For example, should one adopt fractional components for only a part of the overall model (as we did in the examples), simultaneously considers fractional pharmacokinetics, link, active drug and pharmacodynamics? A similar situation, an embarrassment of riches, was first encountered when direct and indirect action models were introduced, and the need for model discrimination arose as a consequence of the novel, larger, class of PKPD models [43]. In addition, there are of course a much larger number of PD models used in the literature [34] than the basic classes we considered here. Tolerance and feedback models might be in particular potentially useful candidates for the application of fractional dynamics.

In conclusion, in this paper we introduce new families of PKPD models based on the application of fractional calculus to PKPD models. The aim of the paper is not to claim the superiority of fractional dynamics models is respect to standard ones, but it is simply to define the new families, and provide some insights into their qualitative behaviour. The main purpose is to suggest possible applications and stimulate further research, which might, or might not, demonstrate the applicability and importance of PKPD models based on fractional calculus.


This work was supported in part by NIH grants R01 AI50587, GM26696.

Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.

Appendix 1

We report the algorithm used to solve the FDE in the examples reported above. The algorithm (FSDEs) is implemented in R [32], and follows the derivation reported in [30, 31]. Both the FI and FDE algorithms are available from the author upon request.

Input variables
  • f = a p-valued function of p + 1 real variables that defines the right-hand side of the differential equation (10).
  • β (beta) = a vector of parameters used in the function f.
  • α (alpha) = the order of the system of FDE (a positive real number).
  • tol = the precision required in the solution.
  • y0 = an array of dimension equation M99, containing the initial values equation M100equation M101, equation M102.
  • xmax = the upper bound of the interval where the solution is to be approximated (a positive real number).
  • n = the number of time steps that the algorithm is to take (a positive integer).

Main internal variables
  • a, b = arrays of n + 1 real numbers that contain the weights of the corrector and predictor formulae, respectively.

Output variables
  • xc = an array of dimension n that contains the times at which the approximate solutions to Eq. 10 is computed.
  • yc = an array of dimension equation M103 real numbers that contains the approximate solutions to Eq. 10, equation M104
    Table thumbnail
    Body of the Algorithm

Appendix 2

We discuss some the properties of the numerical approximation to the Riemann–Liouville integral with particular regard to its relationship with the active drug component the PD models we propose, see Eq. 15 and Table 1.

Figure 6 shows a visualization of the quadrature weights used for fractional integration (similar weights are used for the predictor in the FDE). The line segments displayed represent scaled weights (to retain an easy interpretation of the results), averaged over each subinterval of integration, for an arbitrary integration interval. For β = 1, i.e., equation M105, the average weight is constant, showing how the algorithm reduces to classic trapezoidal integration. Previous drug concentrations are simply “summed” to obtain the exposure to drug. However, as pointed out in [26] there is a quite striking difference for the cases 0 < β < 1 and 1 < β < 2.

Fig. 6
Scaled weights of quadrature for the numerical approximation of the Riemann–Liouville fractional integral (Eq. 5) and active drug equation M106)

When 0 < β < 1, the weights of earlier values contribute less to the overall solution than the more recent ones, that is fractional integration exhibits a long-term memory loss in respect to standard integration, and the smaller the value of β the greater the degree of long-term memory loss. In terms of PD models 0 < β < 1 corresponds to hypo-exposure models, and the weights indicate how recent drug concentrations contribute more to the value of “active drug” than previous ones; to the limit when β = 0, i.e., equation M107, the effect only depends on current drug concentration, so that for a direct action model we are back to the standard memory-less relationship between effect and concentration in the effect compartment.

In contrast, when 1 < β < 2, the earlier values of the integrand will contribute more to the overall solution than will the more recent ones. Fractional integration therefore exhibits a short-term memory loss, in respect to standard integration. Furthermore, the greater the value of β the more pronounced the short-term memory loss will be. The case 1 < β < 2, corresponds to hyper-exposure PD models, providing a representation in which the most recent drug concentrations are less important than the concentrations at the beginning of the administration to determine the effect.


1A system of m fractional differential equations is said to be commensurate when each FDE is of the same fractional order, α; non-commensurate when the fractional orders, α1, α2, …, αm, are distinct for each equation. For the PD models we propose α, β, and δ can be distinct.


1. Caputo M. Linear models of dissipation whose q is almost frequency independent—II. Geophys J Roy Astron Soc. 1967;13:529–539.
2. Caputo M, Mainardi F. Linear models of dissipation in anelastic solids. Rivista del Nuovo Cimento. 1971;1:161–198. doi: 10.1007/BF02820620. [Cross Ref]
3. Sokolov IM, Klafter J, Blumen A. Fractional kinetics. Physics Today. 2002;55:48–54. doi: 10.1063/1.1535007. [Cross Ref]
4. Oldham KB, Spanier J. The fractional calculus. San Diego, CA: Academic Press; 1974.
5. II RJM, Hall MW. Differintegral interpolation from bandlimited signals samples. IEEE Trans Acoust Speech Signal Process. 1981;29:872–877. doi: 10.1109/TASSP.1981.1163636. [Cross Ref]
6. Olmstead WE, Handelsman RA. Diffusion in a semi-infinite region with nonlinear surface dissipation. SIAM Rev. 1976;18:275–291. doi: 10.1137/1018044. [Cross Ref]
7. Chern JT (1993) Finite element modeling of viscoelastic materials on the theory of fractional calculus. Ph.D. thesis, University Microfilms No. 9414260
8. Gaul L, Klein P, Kempfle S. Damping description involving fractional operators. Mech Syst Signal Process. 1991;5:81–88. doi: 10.1016/0888-3270(91)90016-X. [Cross Ref]
9. Glockle WG, Nonnenmacher TF. A fractional calculus approach to self-similar protein dynamics. Biophys J. 1995;68:46–53. doi: 10.1016/S0006-3495(95)80157-8. [PubMed] [Cross Ref]
10. Metzler R, Schick W, Kilian H-G, Nonnenmacher TF. Relaxation in filled polymers: a fractional calculus approach. J Chem Phys. 1995;103:7180–7186. doi: 10.1063/1.470346. [Cross Ref]
11. Bagley RL, Torvik PJ. A theoretical basis for the application of fractional calculus to viscoelasticity. J Rheol. 1983;27:201–210. doi: 10.1122/1.549724. [Cross Ref]
12. Unser M, Blu T. Fractional splines and wavelets. SIAM Rev. 2000;42:43–67. doi: 10.1137/S0036144598349435. [Cross Ref]
13. Forster B, Blu T, Unser M. Complex b-splines. Appl Comp Harm Anal. 2006;20:261–282. doi: 10.1016/j.acha.2005.07.003. [Cross Ref]
14. Podlubny I (1994) Fractional-order systems and fractional-order controllers. Tech. Report UEF-03-94, Institute for Experimental Physics, Slovak Academy of Sciences
15. Xin C, Fawang L. Numerical simulation of the fractional-order control system. J Appl Math Comput. 2007;23:229–241. doi: 10.1007/BF02831971. [Cross Ref]
16. El-Sayed AMA, Rida SZ, Arafa AAM. On the solutions of time-fractional bacterial chemotaxis in a diffusion gradient chamber. Int J Nonlin Sci. 2009;7:485–492.
17. Nonnenmacher TF, Metzler R. On the Riemann–Liouville fractional calculus and some recent applications. Fractals. 1995;3:557–566. doi: 10.1142/S0218348X95000497. [Cross Ref]
18. Podlubny I. Geometric and physical interpretation of fractional integration and fractional differentiation. Fract Calc Appl Anal. 2002;5:367–386.
19. Dokoumetzidis A, Macheras P. Fractional kinetics in drug absorption and disposition processes. J Pharmacokinet Pharmacodyn. 2009;36:165–178. doi: 10.1007/s10928-009-9116-x. [PubMed] [Cross Ref]
20. Popovic JK, Atanackovic MT, Pilipovic AS, Rapaic MR, Pilipovic S. A new approach to the compartmental analysis in pharmacokinetics: fractional time evolution of diclofenac. J Phamacokinet Pharmacodyn. 2010;37:119–134. doi: 10.1007/s10928-009-9147-3. [PubMed] [Cross Ref]
21. Verotta D. Fractional compartmental models and multi-term mittag-leffler response functions. J Phamacokinet Pharmacodyn. 2010;37:209–215. doi: 10.1007/s10928-010-9155-3. [PMC free article] [PubMed] [Cross Ref]
22. Miller KS, Ross B. An introduction to the fractional calculus and fractional differential equations. New York: Wiley; 1993.
23. Podlubny I. Fractional differential equations. San Diego: Academic Press; 1999.
24. Davis PJ. Leonhard Euler’s integral: a historical profile of the gamma function. Am Math Monthly. 1959;66:849–869. doi: 10.2307/2309786. [Cross Ref]
25. Lindstrom MJ, Bates DM. Nonlinear mixed effect models for repeated measures data. Biometrics. 1990;46:673–687. doi: 10.2307/2532087. [PubMed] [Cross Ref]
26. Diethelm K, Ford NJ. Analysis of fractional differential equations. J Math Anal Appl. 2002;265:229–248. doi: 10.1006/jmaa.2000.7194. [Cross Ref]
27. Kilbas AA, Trujillo JJ. Differential equations of fractional order: methods, results and problems I. Appl Anal. 2001;78:153–192. doi: 10.1080/00036810108840931. [Cross Ref]
28. Magin RL. Fractional calculus in bioengineering—part 2. Crit Rev Biomed Eng. 2004;32:105–193. doi: 10.1615/CritRevBiomedEng.v32.i2.10. [PubMed] [Cross Ref]
29. Magin RL. Fractional calculus in bioengineering—part 3. Crit Rev Biomed Eng. 2004;32:195–377. doi: 10.1615/CritRevBiomedEng.v32.i34.10. [PubMed] [Cross Ref]
30. Diethelm K, Ford NJ, Freed AD, Luchko Y. Algorithms for the fractional calculus: a selection of numerical methods. Comput Methods Appl Mech Eng. 2005;194:743–773. doi: 10.1016/j.cma.2004.06.006. [Cross Ref]
31. Diethelm K, Ford JM, Ford NJ, Weilbeer M. Pitfalls in fast numerical solvers for fractional differential equations. J Comput Appl Math. 2006;186:482–503. doi: 10.1016/ [Cross Ref]
32. The R project for statistical computing. (2006) [PubMed]
33. Gorenflo R. Fractional calculus: some numerical methods. Wien: Springer; 1997.
34. Csajka C, Verotta D. Pharmacokinetics and pharmacodynamics modeling: history and perspectives. J Pharmacokin Pharmacodyn. 2006;33:227–279. doi: 10.1007/s10928-005-9002-0. [PubMed] [Cross Ref]
35. Segre G. Kinetics of interaction between drugs and biological systems. Farmaco Sci. 1968;23:907–918. [PubMed]
36. Dayneka NL, Garg V, Jusko WJ. Comparison of four basic models of indirect pharmacodynamic responses. J Pharmacokin Biopharm. 1993;21:457–478. doi: 10.1007/BF01061691. [PMC free article] [PubMed] [Cross Ref]
37. Erdélyi A, Magnus W, Oberhettinger F, Tricomi FG. §8.2.—Bateman Manuscript Project: higher transcendental functions, vol 2. New York: McGraw-Hill; 1955.
38. Pitsiu M, Gries J, Benowitz N, Gourlay SG, Verotta D. Modeling nicotine arterial-venous differences to predict arterial concentrations and input based on venous measurements: application to smokeless tobacco and nicotine gum. J Pharmacokin Pharmacodyn. 2002;29:383–402. doi: 10.1023/A:1020957208071. [PubMed] [Cross Ref]
39. Zamek-Gliszczynski MJ, Kalvass JC, Pollack GM, Brouwer KLR. Relationship between drug/metabolite exposure and impairment of excretory transport function. Drug Metab Dispos. 2009;37:386–390. doi: 10.1124/dmd.108.023648. [PubMed] [Cross Ref]
40. Boeckmann AJ, Beal SL, Sheiner LB (2006) NONMEM VI user’s guides technical report. Project Group, University of California at San Francisco, San Francisco
42. Jacquez JJ. Parameter estimation: local identifiability of parameters. Am J Physiol. 1990;258:E727–E736. [PubMed]
43. Hutmacher MM, Mukherjee D, Kowalski KG, Jordan DC. Collapsing mechanistic models, maintaining their interpretations, and the implications on substantive inference: an application to dose selection for proof of concept of a selective irreversible antagonist. J Pharmacokin Pharmacodyn. 2005;32:501–520. doi: 10.1007/s10928-005-0052-0. [PubMed] [Cross Ref]

Articles from Springer Open Choice are provided here courtesy of Springer