|Home | About | Journals | Submit | Contact Us | Français|
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.
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  and Caputo and Mainardi  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 . Other examples are of applications are in diffusion processes , signal processing , diffusion problems . More recent applications are mainly in physics: finite element implementation of viscoelastic models , mechanical systems subject to damping , relaxation and reaction kinetics of polymers , relaxation in filled polymer networks , viscoelastic materials , although there are recent applications in splines and wavelets [12, 13], control theory [14, 15], and biology  (bacterial chemotaxis). Surveys with collections of applications can also be found in Nonnenmacher and Metzler , and Podlubny . The suggestion of possible PK applications is reported in [19, 20], and generalized in , 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.
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 , 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 , here we give a brief overview.
Although mathematical modelers dealing with dynamical systems are very familiar with derivatives of integer order, and their inverse operation, integrations, they are generally much less so with fractional-order derivatives, for example One way to formally introduce fractional derivatives proceeds from the repeated differentiation of an integer power p:
For p a real number, repeated differentiation gives
where the gamma functions  replace the factorials. The gamma function allows for a generalization to a non-integer order of differentiation, with α real,
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
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
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 is defined by
where m is the smallest integer greater than α, ( integer) is the classical differential operator, and f(t) is required to be continuous and α-times differentiable in t. The operator is named after Caputo , 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:
which is important in PK and PD applications to, e.g., preserve linearity in respect to inputs; that it commutes:
which is important to compute partial derivatives in respect to parameters for model fitting, or for future applications involving mixed effect models ; and that, as the standard differential operator, it possesses the desirable property that:
for any constant c, e.g. the Caputo derivative for steady-state values is zero, as it should for physical realism.
Having defined 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
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
where 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 .
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 , 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 . 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.
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 . 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, and Caputo differential equations, are the operators which are used to generalize the standard PKPD models, as described below.
A pharmacodynamic model is traditionally divided in two parts: a link model and a transduction model (see  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 ). 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) . See  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.
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:
where represents the drug input into the central compartment and is the elimination rate constant from the effect compartment. In general the solution of Eq. 12 takes the form of the convolution:
where and is the Mittag–Leffler function :
The Mittag–Leffler function is related to the monoexponential by the relationship and it represents the unit-response function for fractional kinetics (see  for general expressions of the unit-response functions for arbitrary fractional compartmental structures).
Note that when Ce0 = 0, as for the standard effect compartment model, but it also holds a fact that might induce a-posteriori identifiability problems.
In the following example 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 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.
In the standard direct action model the effect at time t, is expressed by an arbitrary (memory-less) function of drug concentration in the effect site at time 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:
For a direct action model the effect is simply related to A(t) by a non-linear memory-less function G:
Equation 16 for β = 0 represents a standard direct action model, the acting drug is drug concentration in the effect site: for β = 1 it represents an effect related to the exposure,, where ; for 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 we have the cases of hyper-exposure dependencies.
The upper panel of Fig. 2 shows 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 for , while for β = 1 .
The lower panel of Fig. 2 shows the corresponding effect profiles when:
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, , but the return is progressively delayed. For β = 1, the effect, reaches the maximum , and does not return to baseline.
The upper panel of Fig. 3 shows 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 is, for β > 1 we have. The lower panel of Fig. 3, shows the effect profiles computed according to (17), where now, for β > 1, .
Dyneka et al. proposed four indirect response models  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:
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
where drug stimulates or inhibits the elimination rate of Y(t). Fractional dynamics can be immediately introduced by substituting , or , 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:
when kin = 1, kout = 1, , Ce(t) is as in the previous examples. The upper panel of Fig. 4 shows the effects profiles corresponding to model (20), for 1. The main features are the quicker onset and slower decay when , dashed lines, in respect to the standard model , solid line. The lower panel of Fig. 4 shows effects profiles for . Now when the onset of the effect is delayed, and damped oscillations are introduced on the way to the return to baseline kin/kout = 1.
Using a similar layout, Fig. 5 shows the effect profiles corresponding to an inhibition of elimination rate:
where all the parameters are as in Fig. 4, but 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.
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 ,  extends the approach to commensurate1 two compartmental models, while  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 , 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.
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.
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 .
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.)
The use of fractional dynamics as in models (20) and (21), allows to smoothly transition from models with exponential or sub-exponential tails (), 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, , 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.
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 . 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 .
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  or NLME  for model identification. In particular the algorithm used here  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 , 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 ) 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 . In addition, there are of course a much larger number of PD models used in the literature  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.
We report the algorithm used to solve the FDE in the examples reported above. The algorithm (FSDEs) is implemented in R , and follows the derivation reported in [30, 31]. Both the FI and FDE algorithms are available from the author upon request.
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., , 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  there is a quite striking difference for the cases 0 < β < 1 and 1 < β < 2.
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., , 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.