recent papers by   and   show the application of fractional calculus to pharmacokinetics ( pk ) .
leffler function : \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ e_{\alpha } \left ( z \right ) = \sum\limits_{i = 0}^{\infty } { { \frac{{z^{i } } } { { \upgamma \left ( { \alpha i + 1 } \right ) } } } } $ $ \end{document}in place of the mono - exponential:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ e^{z } = \sum\limits_{i = 0}^{\infty } { { \frac{{z^{i } } } { i ! } } } = \sum\limits_{i = 0}^{\infty } { { \frac{{z^{i } } } { { \upgamma \left ( { i + 1 } \right ) } } } } = e_{1 } \left ( z \right ) $ $ \end{document}if drug is given as a bolus dose in a venous site , drug concentration in plasma at a time t after the dose administration , c(t ) , can be represented using a response function of the form:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ c(t ) = dose \ , \left [ { \theta e_{\alpha } \left ( { - \lambda t^{\alpha } } \right ) } \right ] $ $ \end{document}where dose is the amount of drug given ,  and 1 have the interpretation of the elimination rate constant and the reciprocal of volume of distribution of the plasma compartment , respectively , if model ( 3 ) is seen as solution corresponding to a single compartment model described by a ( caputo ) fractional differential equation of order  , see , e.g. , [ 5 , 6 ] .
the paper of   shows solutions for some specific two- and three - compartmental structures described by fractional order kinetics .
however , as pointed out by , while the connection between response function and compartmental structure is immediate for the single compartment case , this is less so for the case of multi - compartmental ones .
the purpose of this communication is to discuss this connection , and to do so we will ( 1 ) clarify the distinction between different types of systems of fractional differential equations , in particular discussing the difference between commensurable and non - commensurable ones , ( 2 ) show solutions for the corresponding response functions , and ( 3 ) discuss their application to the modeling of pk data .
commensurate fractional order linear systems are described by a system of linear fractional differential equations ( fde ) of the form : \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ { \frac{{d_{\blacktriangle } ^{{^{\alpha } } } } } { { { \text{d}}t}}}{\mathbf{x}}(t ) = \left ( { \begin{array}{*{20}c } { { \frac{{d_{\blacktriangle } ^{{^{\alpha } } } } } { { { \text{d}}t}}}x_{1 } ( t ) } \\ \ldots \\ { { \frac{{d_{\blacktriangle } ^{{^{\alpha } } } } } { { { \text{d}}t}}}x_{m } ( t ) } \\ \end{array } } \right ) = \left ( { \begin{array}{*{20}c } { a_{11 } } & \ldots & { a_{1 m } } \\ \ldots & \ldots & \ldots \\ { a_{m1 } } & \ldots & { a_{mm } } \\ \end{array } } \right){\mathbf{x}}(t ) + { \mathbf{f}}(t ) = a{\mathbf{x}}(t ) + { \mathbf{f}}(t ) $ $ \end{document}with initial conditions\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ { \mathbf{x}}(0 ) = { \mathbf{x}}_{0 } $ $ \end{document } , where \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ d_{\blacktriangle } ^{{^{\alpha } } } $ $ \end{document } is the caputo fractional differential operator   of order   >  0 , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ { \mathbf{f}}(t ) $ $ \end{document}is the ( vector valued ) input function to the system .
these systems are called commensurate because all the differential equations are of the same fractional order , . as a consequence t can be shown that the solution to the system of fde   represents the entire state of the system at any given time .
in particular , a compartmental system can be obtained , for 0  <      1 , exactly as for a system of ode , by introducing mass - balance constrains:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ a_{ji } \ge 0,\quad i \ne j $ $ \end{document}and the following:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$\begin{array}{*{20}c } { \text{a}}_{ii } \le 0   \\ \left| { { \text{a}}_{ii } } \right|\ge \sum\limits_{\begin{array}{l } j = 1 \\ j \ne i\end{array } } ^{m } { a_{ji } } \hfill \\ \end{array } $ $ \end{document}which guarantee that all states are non - negative . it can be shown ( see e.g. ) that the solutions to the system of linear fde ( 4 ) depend on the eigenvalues of its characteristic equation , that is , in the laplace domain , \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \det \left ( { b(s ) = a(s ) - si } \right ) $ $ \end{document}. in particular if the eigenvalues are real and distinct the solution to eq . 
4 takes the form1:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ { \mathbf{x}}(t ) = b_{1 } { \mathbf{u}}_{1}^{(1 ) } e_{\alpha } ( \lambda_{1 } t^{\alpha } ) + b_{2 } { \mathbf{u}}_{2}^{(2 ) } e_{\alpha } ( \lambda_{1 } t^{\alpha } ) + \cdots + b_{m } { \mathbf{u}}_{m}^{(m ) } e_{\alpha } ( \lambda_{1 } t^{\alpha } ) $ $ \end{document}where b1 , b2 ,  , bm are constants , 1 , 2 ,  , m and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ { \mathbf{u}}_{1}^{(1 ) } , { \mathbf{u}}_{2}^{(2 ) } , \ldots , { \mathbf{u}}_{m}^{(m ) } $ $ \end{document } are the eigenvalues and eigenvectors of the characteristic equation for ( 4 ) .
it is immediate from ( 7 ) that for a bolus input in the j - th compartment the solution for drug concentration in the same compartment takes the form:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ c_{j } ( t ) = dose\left [ { \theta_{1 } e_{\alpha } ( \lambda_{1 } t^{\alpha } ) + \theta_{2 } e_{\alpha } ( \lambda_{2 } t^{\alpha } ) + \cdots + \theta_{m } e_{\alpha } ( \lambda_{m } t^{\alpha } ) } \right ] $ $ \end{document}which establishes a direct connection with the familiar multi - exponential response function corresponding to ordinary multi - compartment linear systems with distinct eigenvalues .
a non - commensurate fractional order linear system is described by : \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ \left ( { \begin{array}{*{20}c } { { \frac{{d_{\blacktriangle } ^{{^{{\alpha_{1 } } } } } } } { { { \text{d}}t}}}x_{1 } ( t ) } \\ \ldots \\ { { \frac{{d_{\blacktriangle } ^{{^{{\alpha_{m } } } } } } } { { { \text{d}}t}}}x_{m } ( t ) } \\ \end{array } } \right ) = a{\mathbf{x}}(t ) + { \mathbf{f}}(t ) $ $ \end{document}where now 1 , 2 ,  , m are distinct ( real positive ) numbers indicating the fractional order for each equation . as remarked by , in reference to the systems of compartments shown in ( 2 ) , the equations in   do not satisfy mass - balance even if conditions ( 5 ) are satisfied , and in general the solution to the system of fde ( 9 ) does not represent the states of the system .
non - negativity is also no longer guaranteed by the relationships ( 6 ) ( and there are non - trivial issues associated with demonstrating the stability , observability and reach - ability of such systems , see . ) numerical methods must be employed to find the solution to ( 9 ) , since a close form solution equivalent to ( 7 ) does not exist [ 10 , 11 ] .
however , solutions can be obtained if is assumed that the fractional orders are rational numbers , that is i  =  pi / qi where pi , qi are integers , i  =  1 ,  , m ( 12).2 the mathematics necessary to obtain the general solution are quite involved , and for the purpose of this paper we only show a subset of the possible solutions , in particular for c(t ) ( see [ 1113 ] for more general results ) . for a non - commensurate system
, it can be shown that a solution for drug concentration in the j - th compartment takes the form:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ c_{j } ( t ) = dose\left\ { { t^{{\gamma - \alpha_{j } } } \left [ { \theta_{1 } e_{{\gamma , \gamma - \alpha_{j } + 1 } } ( \lambda_{1 } t^{\gamma } ) + \theta_{2 } e_{{\gamma , \gamma - \alpha_{j } + 1 } } ( \lambda_{2 } t^{\gamma } ) + \cdots + \theta_{m } e_{{\gamma , \gamma - \alpha_{j } + 1 } } ( \lambda_{m } t^{\gamma } ) } \right ] } \right\ } $ $ \end{document}where   =  1/q , q  =  m.c.d(q1 ,  , qm ) , and \documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ e_{\alpha , \beta } \left ( z \right ) = \sum\limits_{i = 0}^{\infty } { { \frac{{z^{i } } } { { \upgamma \left ( { \alpha
i + \beta } \right ) } } } } $ $ \end{document } is the two - parameters mittag  leffler function ( 5 ) .
one notices two important facts : first , this solution only depends on the fractional order for compartment j , j , and  , second in direct analogy to the solution ( 8) above , the mittag 
leffler function exponents are determined by the eigenvalues of the characteristic equation for the system .
we now have the ingredients to show a simple example of applications of multi - terms mittag  leffler response functions to fit pk data .
we consider the case m  =  2 , which for a standard ode system generates the response function:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ c(t ) = dose\left ( { \theta_{1 } e^{{\lambda_{1 } t } } + \theta_{2 } e^{{\lambda_{2 } t } } } \right ) $ $ \end{document}for a commensurate fde system obtains:\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ c(t ) = dose\left ( { \theta_{1 } e_{\alpha } ( \lambda_{1 } t^{\alpha } ) + \theta_{2 } e_{\alpha } ( \lambda_{2 } t^{\alpha } ) } \right ) $ $ \end{document}and for a non - commensurate fde system\documentclass[12pt]{minimal } 
				 \usepackage{amsmath } 
				 \usepackage{wasysym } 
				 \usepackage{amsfonts } 
				 \usepackage{amssymb } 
				 \usepackage{amsbsy } 
				 \usepackage{mathrsfs } 
				 \usepackage{upgreek } 
				 \setlength{\oddsidemargin}{-69pt } 
				 \begin{document}$$ c(t ) = dose \ , t^{\gamma - \alpha } \left ( { \theta_{1 } e_{\gamma , \gamma - \alpha + 1 } ( \lambda_{1 } t^{\gamma } ) + \theta_{2 } e_{\gamma , \gamma - \alpha + 1 } ( \lambda_{2 } t^{\gamma } ) } \right ) $ $ \end{document}the parameters 1 , 2 , 1 , 2 ,  ,  are estimated from the data , with the constraints 1 , 2  >  0 , 1 , 2  <  0 , and 0  <   ,     1 , which guarantee that eq .  13 is non - negative and non - increasing ( strictly monotone ) for t    0.3 to evaluate the single and two - parameters mittag  leffler function we implemented a fortran 90 version the algorithm reported in .
we used the computer program nonmem   to obtain the parameters estimates . as a check , fig . 
11 of ,   =  1 and 0  <      2 .
contrary to  , which has a strong influence on the overall shape of the curve , the parameter  has its most pronounced influence on the value of the function at t  =  0.fig . 
1the mittag  leffler function e,(t ) for   =  1 and     ( 0 , 2 ) the mittag 
leffler function e,(t ) for   =  1 and     ( 0 , 2 ) figure  2 shows the fit of models ( 11 ) ( solid line ) , ( 12 ) ( widely dashed line ) , and ( 13 ) ( dashed line ) , to error corrupted data simulated using an eight compartments mammillary model .
note the added flexibility introduced by use of a sum of single- and two - parameters mittag  leffler functions in respect to exponentials : the values of minus twice log - likelihood for the fit of the simulated data were 316.895 , 381.974 , and 406.354 , for models ( 11)(13 ) , respectively .
of course , this is just an example to show the feasibility of the approach : for this simulation , a sum of exponentials would fit the simulated data perfectly well.fig .  2simulated data ( circles ) with superimposed the fit of the response function for a second order ordinary system , eq .  11 in the text ( solid line ) , commensurate fde , eq .  12 ( widely dashed line ) , and non - commensurate fde , eq
.  13 ( dashed line ) simulated data ( circles ) with superimposed the fit of the response function for a second order ordinary system , eq .  11 in the text ( solid line ) , commensurate fde , eq .  12 ( widely dashed line ) , and non - commensurate fde , eq .  13 ( dashed line )
following up on the papers of   and , and the commentary by , the first purpose of this commentary is to further clarify the nature of systems of fde , and in particular to point out the distinction between commensurate and non - commensurate ones .
commensurate systems of fde have a direct relationship with system of ode , and in particular when formulated in terms of compartmental models ( that is , satisfying mass balance and non - negativity constraints , see eqs .  5 ,
non - commensurate systems of fde do not , in general , represent the state of a system . leaving to the side
the issue of what exactly system ( 9 ) represents , one is still justified in using it as a black - box type model for single - input / single - output experiments , as long as physical constraints ( non - negativity in particular ) are satisfied .
to do so response functions are a convenient tool , and we show that , for commensurate and non - commensurate fde , relatively simple ones can be derived which satisfy such requirements .
the solutions for a system of commensurate fde takes the form of the sum of mittag  leffler functions , with a single parameter  , eq . 
8 , while solutions for a system of non - commensurate fde can be expressed by a sum of two - parameters mittag  leffler functions , such as eq . 
, one can establish a direct analogy between the familiar sum of exponentials used in pk , and importantly all the relationships between compartmental transfer rate constants and the intercepts and exponents of the corresponding response functions found in classic textbooks on pk [ 17 , 18 ] , carry forward to fractional differential equations . in conclusion , while insight into the physiological interpretability of fde system might be gained in the future , and the formulation of non - commensurate systems of fde to represent the states of a system ( required by , e.g. , a physiological flow model ) requires further investigation , the response functions ( 8) and ( 10 ) can be used to investigate the existence of pk data sets which might actually show complex fractional kinetics.4
 systems of fractional differential equations ( sfde ) have been increasingly used to represent physical and control system , and have been recently proposed for use in pharmacokinetics ( pk ) by ( j pharmacokinet pharmacodyn 36:165178 , 2009 ) and ( j phamacokinet pharmacodyn , 2010 ) . 
 we contribute to the development of a theory for the use of sfde in pk by , first , further clarifying the nature of systems of fde , and in particular point out the distinction and properties of commensurate versus non - commensurate ones . 
 the second purpose is to show that for both types of systems , relatively simple response functions can be derived which satisfy the requirements to represent single - input / single - output pk experiments . 
 the response functions are composed of sums of single- ( for commensurate ) or two - parameters ( for non - commensurate ) mittag  leffler functions , and establish a direct correspondence with the familiar sums of exponentials used in pk . 