Modeling and administration scheduling of fractional-order pharmacokinetic systems

TL;DR

This paper explores fractional-order models for pharmacokinetics, compares simulation methods, and proposes an optimal drug administration schedule, demonstrating its effectiveness on patient data.

Contribution

It introduces simulation techniques for fractional-order pharmacokinetic models and develops an optimal drug scheduling scheme evaluated on patient data.

Findings

Simulation methods vary in accuracy and complexity

Optimal scheduling improves drug administration outcomes

Modeling captures nonlinear pharmacokinetic phenomena

Abstract

Fractional-order dynamical systems were recently introduced in the field of pharmacokinetics where they proved powerful tools for modeling the absorption, disposition, distribution and excretion of drugs which are liable to anomalous diffusion, deep tissue trapping and other nonlinear phenomena. In this paper we present several ways to simulate such fractional-order pharmacokinetic models and we evaluate their accuracy and complexity on a fractional-order pharmacokinetic model of Amiodarone, an anti-arrhythmic drug. We then propose an optimal administration scheduling scheme and evaluate it on a population of patients.

Tables

Table 1. Table 1: Errors using the Padé approximation.

Order	$\Vert e_1\Vert$	$\Vert e_2\Vert$	${\Vert e_1\Vert }_{\infty }$	${\Vert e_2\Vert }_{\infty }$
$\left[2/3\right]$	$2.833\cdot {10}^{-4}$	$1.907\cdot {10}^{-4}$	$0.0015$	$0.001$
$\left[3/4\right]$	$1.105\cdot {10}^{-4}$	$0.0059$	$6.094\cdot {10}^{-4}$	$0.0113$
$\left[4/5\right]$	$4.514\cdot {10}^{-5}$	$2.406\cdot {10}^{-5}$	$3.076\cdot {10}^{-4}$	$1.774\cdot {10}^{-4}$
$\left[5/6\right]$	$2.327\cdot {10}^{-5}$	$2.685\cdot {10}^{-5}$	$1.752\cdot {10}^{-4}$	$1.976\cdot {10}^{-4}$

Table 2. Table 2: Errors using the Oustaloup approximation.

${\omega }_b$	${\omega }_h$	$N$	$\Vert e_1\Vert$	$\Vert e_2\Vert$	${\Vert e_1\Vert }_{\infty }$	${\Vert e_2\Vert }_{\infty }$
${10}^{-2}$	${10}^3$	$8$	$0.0023$	$0.0043$	$0.0228$	$0.0025$
${10}^{-2}$	${10}^4$	$20$	$5.744\cdot {10}^{-4}$	$0.0045$	$0.0054$	$0.0027$
${10}^{-3}$	${10}^3$	$8$	$0.0033$	$5.084\cdot {10}^{-4}$	$0.0332$	$6.555\cdot {10}^{-4}$
${10}^{-3}$	${10}^4$	$20$	$7.451\cdot {10}^{-4}$	$4.597\cdot {10}^{-4}$	$0.0074$	$2.765\cdot {10}^{-4}$

Table 3. Table 3: Errors using the Matsuda-Fujii approximation with s k = β α k subscript 𝑠 𝑘 superscript 𝛽 subscript 𝛼 𝑘 s_{k}=\beta^{\alpha_{k}} .

$\beta$	${\alpha }_k$	$\Vert e_1\Vert$	$\Vert e_2\Vert$	${\Vert e_1\Vert }_{\infty }$	${\Vert e_2\Vert }_{\infty }$
$2$	$-1:10$	$7.01\cdot {10}^{-5}$	$3.162\cdot {10}^{-4}$	$4.82\cdot {10}^{-5}$	$2.111\cdot {10}^{-4}$
$2$	$1:10$	$0.0016$	$0.0036$	$0.0009$	$0.0032$
$2.3$	$-1:11$	$0.0002$	$0.0032$	$0.0018$	$0.003$
$3$	$1:10$	$0.0034$	$0.0127$	$0.0017$	$0.0075$

Table 4. Table 4: Approximation errors of fde12

$h$	$\Vert e_1\Vert$	$\Vert e_2\Vert$	${\Vert e_1\Vert }_{\infty }$	${\Vert e_2\Vert }_{\infty }$
${10}^{-2}$	$0.0333$	$0.0223$	$0.0584$	$0.0368$
${10}^{-3}$	$0.0020$	$0.0013$	$0.0083$	$0.0055$
${10}^{-4}$	$1.782\cdot {10}^{-4}$	$1.157\cdot {10}^{-4}$	$7.982\cdot {10}^{-4}$	$5.204\cdot {10}^{-4}$
${10}^{-5}$	$1.669\cdot {10}^{-5}$	$1.113\cdot {10}^{-5}$	$7.412\cdot {10}^{-4}$	$4.824\cdot {10}^{-4}$

Table 5. Table 5: Approximation errors of the Grünwald-Letnikov method

$h$	$h\nu$	$\Vert e_1\Vert$	$\Vert e_2\Vert$	${\Vert e_1\Vert }_{\infty }$	${\Vert e_2\Vert }_{\infty }$
${10}^{-2}$	$3$	$8.165\cdot {10}^{-4}$	$0.0022$	$8.223\cdot {10}^{-4}$	$0.0021$
${10}^{-2}$	$5$	$5.302\cdot {10}^{-4}$	$4.571\cdot {10}^{-4}$	$8.223\cdot {10}^{-4}$	$5.911\cdot {10}^{-4}$
${10}^{-2}$	$7$	$5.126\cdot {10}^{-4}$	$3.640\cdot {10}^{-4}$	$8.223\cdot {10}^{-4}$	$5.911\cdot {10}^{-4}$
${10}^{-3}$	$3$	$6.390\cdot {10}^{-4}$	$0.0022$	$4.735\cdot {10}^{-4}$	$0.0021$
${10}^{-3}$	$5$	$1.454\cdot {10}^{-4}$	$2.760\cdot {10}^{-4}$	$1.727\cdot {10}^{-4}$	$4.057\cdot {10}^{-4}$
${10}^{-3}$	$7$	$5.060\cdot {10}^{-5}$	$3.594\cdot {10}^{-5}$	$8.067\cdot {10}^{-4}$	$5.797\cdot {10}^{-4}$

Equations

$$\frac{\mathrm{d}{A}_1}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{A}_2}{\mathrm{d}t}$$

$$F(x(t),D^{{\alpha }_1}x(t),\dots ,D^{{\alpha }_p}x(t))=0,$$

$$({{}_{\mathrm{rl}}I}^{\alpha }f)(t)=\frac{1}{\Gamma (\alpha )}{\int }_0^t(t-\tau {)}^{\alpha -1}f(\tau )\mathrm{d}\tau ,t\ge 0.$$

$$({{}_{\mathrm{c}}D}^{\alpha }f)(t)={{}_{\mathrm{rl}}I}^{m-\alpha }\frac{{\mathrm{d}}^{m}f(t)}{\mathrm{d}{t}^m}.$$

$$({{}_{\mathrm{gl}}D}^{\alpha }f)(t)=\underset{h\to 0}{\lim }\frac{1}{h^{\alpha }}\sum _{i=0}^{\infty }(-1{)}^i(\genfrac{}{}{0ex}{}{\alpha }{i})f(t-ih),$$

$$G_1(s)$$

$$G_2(s)$$

$$T(s)=\frac{P(s)}{Q(s)},$$

$$H(s)={\alpha }_0+\frac{s-s_0}{{\alpha }_1+\frac{s-s_1}{{\alpha }_2+\frac{s-s_2}{{\alpha }_3+\dots }}}$$

$$H(s)=s^{\alpha },$$

$$\hat{H}(s)=c_0\prod _{k=-N}^N\frac{s+{\omega }_k}{s+{\omega }_k^{\prime }}$$

$${\omega }_k^{\prime }={\omega }_b{\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\frac{k+N+0.5(1+\alpha )}{2N+1}},$$

$${\omega }_k={\omega }_b{\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{\frac{k+N+0.5(1-\alpha )}{2N+1}},$$

$$c_0={\left(\frac{{\omega }_h}{{\omega }_b}\right)}^{-\frac{r}{2}}\prod _{k=-N}^N\frac{{\omega }_k}{{\omega }_k^{\prime }}.$$

$$({{}_{\mathrm{gl}}\Delta }^{\alpha }x{)}_k=\frac{1}{h^{\alpha }}\sum _{i=0}^{\infty }(-1{)}^i\left(\genfrac{}{}{0ex}{}{\alpha }{i}\right)x_{k-i},$$

$$({{}_{\mathrm{gl}}\Delta }_{\nu }^{\alpha }x{)}_k=\frac{1}{h^{\alpha }}\sum _{i=0}^{\nu }(-1{)}^i\left(\genfrac{}{}{0ex}{}{\alpha }{i}\right)x_{k-i},$$

$${{}_{c}D}^{\gamma }x(t)$$

$$x^{(k)}(0)$$

$${}_c{D}^{\gamma }{x}_0$$

$${}_c{D}^{\gamma }{x}_1$$

$$⋮$$

$${}_c{D}^{\gamma }{x}_{q-1}$$

$${}_c{D}^{\gamma }{x}_q$$

$$⋮$$

$${}_c{D}^{\gamma }{x}_{2q-1}$$

$$x(t)=\sum _{k=0}^{m-1}{x}_{0,k}\frac{t^k}{k!}+({{}_{\mathrm{rl}}I}^{\gamma }f)(t,x(t)),$$

$$x(t_{n+1})$$

$$+\sum _{j=1}^n{a}_{j,n+1}f(t_j,x(t_j)).$$

$$x_P(t_{n+1})=T_{m-1}(t_n)$$

$$({{}_{\mathrm{rl}}I}_h^{\gamma }f)(t)\approxeq h^{\gamma }\sum _{j=0}^n{\omega }_{n-j}f(t_j)+h^{\gamma }\sum _{j=0}^s{w}_{n,j}f(t_j),$$

$$f(t)=\underset{T\to \infty }{\lim }\frac{1}{2\pi i}{\int }_{\sigma -iT}^{\sigma +iT}{\text{e}}^{st}F(s)\text{d}s.$$

$$\Vert e_i\Vert$$

$$\Vert e_i{\Vert }_{\infty }$$

$$t_c^{-1}(x_{k+1}-x_k)=A{x}_k+F{{}_{\mathrm{gl}}\Delta }_{\nu }^{1-a}{x}_k+B{u}_k$$

$$A=\left[\begin{array}{cc}-(k_{12}+k_{10})& 0\\ k_{21}& 0\end{array}\right],F=\left[\begin{array}{cc}0& k_{21}\\ 0& -k_{21}\end{array}\right],B=\left[\begin{array}{c}1\\ 0\end{array}\right].$$

$$\displaystyle x_{k+1}=x_{k}+t_{c}\bigg{(}Ax_{k}+\frac{F}{t_{c}^{1-a}}\sum_{j=0}^{\nu}c_{j}^{1-a}x_{k-j}+Bu_{k}\bigg{)},$$

Full text

††thanks: The work of the last author was supported by the KU Leuven Research Council under BOF/STG-15-043.

, , , , , and

1 Introduction

Pharmacokinetic (PK) models are systems of differential equations which simulate the dynamic response of living organisms in terms of the concentration of a drug or any other substance in different compartments (organs) of the body following its administration to the body. PK models can assist in designing effective and safe administration strategies for individual patients or populations thereof. Among the different types of PK modeling, fractional PK models have attracted the interest of researchers in the field because they can model phenomena like anomalous diffusion, deep tissue trapping and diffusion across fractal manifolds, which traditional PK models fail to describe (Dokoumetzidis and Macheras, 2008, 2011; Dokoumetzidis et al., 2010).

However, simulating such dynamics is not as straightforward as with integer-order systems. Analytical solutions are rarely available and, even then, the evaluation of the solution requires a numerical approximation method (Kaczorek, 2011). The availability of accurate discrete-time approximations of the trajectories of such systems is important not only for simulating but also for the design of open-loop or closed-loop (such as with model predictive control) administration strategies (Sopasakis et al., 2015). Here, we compare several approximation methods (such as ones based on the Grünwald-Letnikov operator and the Oustaloup filter) using a specific PK model taken from the recent literature (Dokoumetzidis et al., 2010) and we discuss their accuracy and fitness for control design.

In this paper we provide a survey of different approaches for modeling fractional-order systems which arise in pharmacokinetics. Our discussion revolves around the case study of Amiodarone, an anti-arrhythmic drug that exhibits fractional-order dynamics. We identify three major classes of numerical algorithms for simulating fractional-order systems in the literature: (i) using rational transfer functions, (ii) time-domain methods and (iii) the numerical inverse Laplace approach. We discuss their merits and limitations and we present a comparative assessment regarding precision of various methods. Our goal, however, is to single out a method which is most suitable for controller design. In the last section we formulate an optimal control problem for administration scheduling to confirm our findings.

2 Fractional Pharmocokinetics: Modeling and Simulation

2.1 Fractional-order pharmacokinetics

Amiodarone is an anti-arrhythmic agent which can be administered either intravenously (i.v.) or orally (Kühlkamp et al., 1999). It is well-known for its highly nonlinear non-exponential dynamics and singular long-term accumulation pattern. Dokoumetzidis et al. (2010) modeled the pharmacokinetic distribution of Amiodarone with a fractional compartmental model following a single intravenous and a single oral dose. The compartmental topology of the model is presented in Figure 1 where it is shown that the diffusion from the tissues to the central compartment is governed by a fractional-order dynamics.

Let $A_{1}$ and $A_{2}$ be the amounts of Amiodarone (in $\mathrm{ng}$) in the plasma and the tissues respectively and $u$ be the administration rate (in $\mathrm{ng}/\mathrm{day}$). We assume that the drug is administered directly into the central (plasma) compartment while the control objective is the concentration of the drug in the tissues attains a prescribed value (set-point). The fractional dynamical model we employ reads as follows:

[TABLE]

with $\alpha\in(0,1)$ and ${{}_{\mathrm{c}}\operatorname{D}}^{1-\alpha}$ is the Caputo fractional derivative which is defined in the following section.

2.2 Fractional-order derivatives

Several fractional-order derivatives have been proposed in the literature the most popular of which are the Riemann-Liouville ${{}_{\mathrm{rl}}\operatorname{D}}^{\alpha}$, the Caputo ${{}_{\mathrm{c}}\operatorname{D}}^{\alpha}$ and the Grünwald-Letnikov ${{}_{\mathrm{gl}}\operatorname{D}}^{\alpha}$ derivatives (Samko et al., 1993). These operators are used to formulate fractional-order differential equations, that is functional equations of the form

[TABLE]

where $\operatorname{D}^{\alpha}$ is a generalized derivative of order $\alpha\geq 0$. Typically, the Caputo derivative is used in this context as the initial conditions are easier to postulate.

The generalized Riemann-Liouville fractional-order integral operator of order $\alpha>0$ is given by

[TABLE]

For $\alpha\in\mathbb{R}$ let us denote by $m=\lceil\alpha\rceil$ the smallest natural number $m$ so that $m\geq\alpha$. The following operator is known as the Caputo derivative of order $\alpha$:

[TABLE]

The Grünwald-Letnikov fractional-order derivative is defined as

[TABLE]

where ${\alpha\choose 0}=1$ and $\binom{\alpha}{i}=\prod_{l=0}^{i-1}\tfrac{\alpha-l}{l+1}$.

The Laplace transform of the Caputo derivative of fractional order $\alpha\in(0,1)$ with zero initial conditions is given as $\mathcal{L}\left[({{}_{\mathrm{c}}\operatorname{D}}^{\alpha})(t)\right]=s^{\alpha}F(s)$, where $F(s)$ is the Laplace transform of function $f(t)$. The transfer function $G_{i}(s)=A_{i}(s)/U(s)$ which associates the administration rate $U(s)=\mathcal{L}u(t)$ to the concentrations $\hat{A}_{i}(s)=\mathcal{L}A_{i}(t)$ are:

[TABLE]

with $\alpha=0.587$, $k_{10}=1.4913\,\mathrm{day^{-1}}$, $k_{12}=2.9522\,\mathrm{day^{-1}}$ and $k_{21}=0.4854\,\mathrm{day^{-\alpha}}$. The concentration of Amiodarone in the two compartments with initial condition $A_{1}(0)=0.1\,\mathrm{ng}$ and $A_{2}(0)=0\,\mathrm{ng}$ is shown in Figure 2.

2.3 Solutions of FDEs

There can be identified four types of solutions for fractional-order differential equations: (i) analytical solutions, (ii) approximations in the $s$-domain using integer-order rational transfer functions and (iii) numerical approximation schemes in the discrete time domain, (iv) the numerical inverse Laplace transformation.

2.3.1 Analytical solutions.

Analytical solutions, when available, involve special functions such as the Mittag-Leffler function $\mathcal{E}_{\alpha,\beta}(t)=\sum_{k=0}^{\infty}{t^{k}}/{\Gamma(\alpha k+\beta)}$ whose evaluation requires in turn some numerical approximation scheme. Typically for the evaluation of the this function we resort to solving an FDE numerically (Garrappa, 2015a).

2.3.2 Transfer function approximations.

Rational approximations aim at approximating the transfer function of a fractional-order system — which involves terms of the form $s^{\alpha}$ — by ordinary transfer functions of the form

[TABLE]

where $P$ and $Q$ are polynomials and the degree of $P$ is no larger than the degree of $Q$.

Padé Approximation: The Padé approximation of order $[m/n]$, $m,n\in\mathbb{N}$, at a point $s_{0}$ is rather popular and leads to rational functions with $\deg P=m$ and $\deg Q=n$ (Silva et al., 2006).

Matsuda-Fujii Method: This method consists in interpolating a function $H(s)$, which is treated as a black box, across a set of logarithmically spaced points (Matsuda and Fujii, 1993). By letting the selected points be $s_{k}$, $k=0,1,2,\dots$, the approximation is written as the continued fractions expansion

[TABLE]

where, $\alpha_{i}=\upsilon_{i}(s_{i})$, $\upsilon_{0}(s)=H(s)$, $\upsilon_{i+1}(s)=\frac{s-s_{i}}{\upsilon_{i}(s)-\alpha_{i}}$

Oustaloup’s method: Oustaloup’s method is based on the approximation of a function of the form:

[TABLE]

with $\alpha>0$ by a rational function

[TABLE]

within a range of frequencies from $\omega_{b}$ to $\omega_{h}$ (Oustaloup et al., 2000). The Oustaloup method offers an approximation at frequencies which are geometrically distributed about the characteristic frequency $\omega_{u}=\sqrt{\omega_{b}\omega_{h}}$ — the geometric mean of $\omega_{b}$ and $\omega_{h}$. The parameters $\omega_{k}$ and $\omega_{k}^{\prime}$ are determined via the design formulas (Petrás, 2011)

[TABLE]

Parameters $\omega_{b}$, $\omega_{h}$ and $N$ are design parameters of the Oustaloup method.

Other methods: There are a few more methods which have been proposed in the literature to approximate fractional-order systems by rational transfer functions such as (Charef et al., 1992; Carlson and Halijak, 1964), as well as data-driven system identification techniques (Gao and Liao, 2012). Nice (Liang et al., 2014).

2.3.3 Time domain approximations.

Several methods have been proposed which attempt to approximate the solution to a fractional-order initial value problem in the time domain.

Grünwald-Letnikov: This is the method of choice in the discrete time domain where ${{{}_{\mathrm{gl}}\operatorname{D}}}^{\alpha}f$ is approximated by its discrete time variant

[TABLE]

which is in turn approximated by a discrete operator with finite memory $\nu$

[TABLE]

which is proven to have bounded error with respect to $({{{}_{\mathrm{gl}}\Delta}}^{\alpha}f)_{k}$ (Sopasakis and Sarimveis, 2017).

Numerical integration methods: Fractional-order initial value problems can be solved with various numerical methods such as the Adams-Bashforth-Moulton predictor-corrector (ABMPC) method (Zayernouri and Matzavinos, 2016) and fractional linear multi-step methods (FLMMs) (Garrappa, 2015b). These methods are only suitable for a system of FDE’s in the form

[TABLE]

where $\gamma$ is a rational, and $m=\lceil\gamma\rceil$.

In order to bring (1) in this form, we need to find a rational approximation of two derivatives, $1-\alpha$ and $1$. If we can find a satisfying rational approximation of $1-a\approxeq p/q$, then first order derivative follows trivially. Now, (1) can be written as

[TABLE]

subject to $x_{0}=A_{1}(0)$, $x_{q}=A_{2}(0)$ and $x_{i}(0)=0$ for $i\neq 1,q$, and $\gamma=1/q$. This system is in fact a linear fractional-order system for which analytical solutions are available (Kaczorek, 2011).

The number of states of system (15) is $2q$, therefore, the rational approximation should aim at a small $q$. In our case $1-\alpha=0.413$ can be written as $413/1000$, but then we would need to simulate a fractional-order system with $2000$ states. Instead $1-\alpha$ can be approximated by $19/46$ with error $0.413-19/46=-4.3478\cdot 10^{-5}$ or $216/523$ with error $0.413-216/523=-1.912\cdot 10^{-6}$. Such approximations can be obtained by means of continued fractions expansions of $1-\alpha$. Yet another reason to choose small $q$ is that small values of $\gamma=1/q$ render the system hard to simulate numerically.

Adams-Bashforth-Moulton predictor-corrector (ABMPC): Methods of the ABMPC type have been generalized to solve fractional-order systems. The basic concept is to evaluate $({{{}_{\mathrm{rl}}\operatorname{I}}}^{\gamma}f)(t,x(t))$ by approximating $f$ with appropriately selected polynomials. Solutions of (14) satisfy the following integral representation

[TABLE]

where the first term on right hand side will be denoted with $T_{m-1}(t).$ The integral on the right hand side of the previous equation can be approximated, using an uniformly spaced grid $t_{n}=nh$, by $\frac{h^{\gamma}}{\gamma(\gamma+1)}\sum_{j=0}^{n+1}a_{j,n+1}f(t_{j})$ for suitable coefficients $a_{j,n+1}$ (Diethelm et al., 2002). The numerical approximation of the solution of (14) is

[TABLE]

Unfortunately, the convergence error of ABMPC when $0<\gamma<1$ is $\mathcal{O}(h^{1+\gamma})$, therefore, a rather small step size $h$ is required to attain a reasonable aprpoximation error. A modification of the basic predictor-corrector method with more favorable computational cost is provided in (Garrappa, 2010) for which the MATLAB implementation fde12 is available.

Lubich’s method: Fractional linear multistep methods (FLMM) (Lubich, 1986) are a generalization of linear multistep methods (LMM) for ordinary differential equations. The key idea is to approximate the Riemann-Liouville fractional-order integral operator (3) with a discrete convolution, called convolution quadrature, as

[TABLE]

for $t_{j}=jh,h>0$ where starting ($w_{n,j}$) and quadrature weights ($\omega_{n}$) are independent of $h$. Surprisingly, the latter weights can be constructed from any linear multistep method for arbitrary fractional order $\gamma$ (Lubich, 1986). Furthermore, FLMM constructed this way will inherit the same convergence rate and at least the same stability properties as the original LMM method (Lubich, 1985).

Here we use the MATLAB implementation flmm2 (Garrappa, 2015c) which is based on (Garrappa, 2015d). However, the method does not perform well for small $\gamma$. In our case study simulations we have found that values smaller than $0.1$ give poor results and often do not converge. However, when we used a more crude approximation of the original system with $\gamma=1/5$, flmm2 method was outperforming fde12 in terms of accuracy and has shown excellent stability properties with respect to bigger step size $h$.

2.3.4 Numerical inverse Laplace.

Several numerical inverse Laplace methods provide an approximation of

[TABLE]

for a given transfer function $F(s)$. Numerical methods can be used to directly evaluate the inversion integral (19) for non-rational transfer functions. One of the most popular methods is to convert the inversion integral into a Fourier transform and then approximate it by a Fourier series via trapezoid rule (de Hoog et al., 1982). These methods are considered to be precise for a broad class of function, but computationally demanding. An implementation of the above method is freely available on line (Hollenbeck, 1998).

A somewhat different approach is taken by Valsa and Brančik (1998), where authors approximate $\e^{st}$, the kernel of the inverse Laplace transformation, by $\tfrac{\e^{st}}{1+e^{-2a}e^{2st}}$ and choose $a$ appropriately so as to achieve an accurate inversion.

In general, numerical inversion methods can achieve high precision, but they are not suitable for control design purposes, especially for optimal control problems. Moreover, different methods are suited for various types of problems. An overview of the most popular inversion methods used in engineering practice is given in (Hassanzadeh and Pooladi-Darvish, 2007).

2.4 Assessment

In order to assess the accuracy of each approximation method presented above we introduce the following error indices

[TABLE]

where $T_{s}$ is a fixed simulation time and $e_{i}(t)$ is the difference between the approximate response of the system $\hat{A}_{i}$ and the one estimated by the inverse Laplace method of Valsa and Brančik (1998) with $a=11$ which is considered to be the most accurate.

The pharmacokinetic profile following a single i.v. bolus dose is shown in Figure 2. In Figures 3, 4 and 5 we show the modeling errors $e_{i}(t)$ and in Tables 1, 2 and 3 we show the corresponding total errors. It seems that adequately high precision can be achieved with these methods. However, approximations in the $s$-domain are not suitable for constrained systems since there is no theoretical bound on the approximation error in the time domain.

The errors of fde12 are presented in Figure 6 and Table 4. Using a step size as small as $h=10^{-5}$ fde12 achieves an approximation error which is uniformly lower than $10^{-4}$. In Figure 7 and Table 4 we show the approximation errors of the Grünwald-Letnikov method. It can be seen that the use of a long history is more important for the attainment of high precision compared to a small step size. Most likely as a result of the low value of $\gamma=0.0217$ in (15), flmm2 failed to produce reasonable approximations. We were in fact only able to produce moderately accurate approximations using $1-\alpha\approxeq 2/5$ where $\gamma=0.2$.

3 Administration scheduling

In this section we address the problem of administration scheduling. Our objective is to devise an administration schedule — a sequence of dosages — so that the concentration of Amiodarone in the tissues is close to a desired value, while the concentration in both compartments never exceeds certain safety limits. We also have to account for a limit on allowed drug dose at each time instant. In doing so, we must assume that we are not able to measure drug concentrations during the treatment. All of these requirements and constraints can be elegantly integrated within the framework of constrained optimal control.

3.1 Optimal control

In this section we describe the optimal control problem formulation. We start by discretizing (1) with a sampling time $t_{c}$ yielding

[TABLE]

where $x_{k}=\left[A_{1}(kt_{c})\,\,A_{2}(kt_{c})\right]^{\prime}$. The left hand side of (22) corresponds to the forward Euler approximation of the first-order derivative, and we shall refer to $t_{c}=10^{-2}\,\mathrm{days}$ as the control sampling time. Matrices $A,F$ and $B$ are

[TABLE]

The discrete-time dynamic equations of the system can now be stated as

[TABLE]

where $c_{j}^{\alpha}=(-1)^{j}\binom{\alpha}{j}$. By augmenting the system with past values as $\tilde{x}_{k}=(x_{k},x_{k-1},\dots,x_{k-\nu+1})$ we can rewrite (24) as a finite-dimension linear system

[TABLE]

Matrices $\hat{A}$ and $\hat{B}$ are straightforward to derive and are given in (Sopasakis and Sarimveis, 2017). The therapeutic session will last for $N_{d}=Nt_{c}=7\,\mathrm{days}$ in total, where $N$ is called the prediction horizon. It is not realistic to administer the drug to the patient too frequently, so we assume that the patient is to receive their treatment every $t_{d}=0.5\,\mathrm{days}$. The administration schedule must ensure that the concentration of drug in all compartments never exceeds the minimum toxic concentration limits while tracking the prescribed reference value as close as possible. To this aim we postulate the following constrained optimal control problem.

[TABLE]

In the above formulation $x_{\mathrm{ref},k}$ is the desired drug concentration at time $k$ and operator ′ denotes vector transposition. Any deviation from set point is penalized by weight matrix $Q=\mathrm{diag}([0\;1])$. Note that we are tracking only the second state. Our underlying GL model has a relative history of $t_{c}\nu=5\,\mathrm{days}$. Optimal drug concentrations are denoted by $u^{\star}_{k}$, for $k=0,\dots,N_{d}-1$ and they correspond to dosages administered intravenously at times $kt_{d}$. In the optimal control formulation we have implicitly assumed that $t_{d}$ is an integer multiple of $t_{c}$, which is not restrictive since $t_{c}$ can be chosen arbitrarily. Finally, we can recognize that problem (26) is a standard quadratic problem that can be readily solved. In our simulations, we have used YALMIP (Löfberg, 2004) to model the problem and MOSEK (ApS, 2015) as the underlying solver to calculate the solution.

3.2 Simulations

To argue for the soundness and the applicability of our approach in real-world scenarios, we will apply the optimal drug dosage schedule to a more precise model than (25). For this purpose we will use fde12 solver. As evident from the results in the previous section, for sufficiently small solver time $h_{sol}$, we can have a realistic simulation of the system.

After solving the optimal control problem we applied the optimal sequence to the FDE simulator fde12. Results are shown in Figures 8 and 9. It can be seen that open loop predictions of the GL model and fde12 simulation show high agreement. This should not be surprising considering that all of the parameters describing the patient are nominal and, additionally, shows that the GL model is of good quality.

Next, we consider multiple patients that are characterized by perturbed parameters. We will simulate the behaviour of $100$ different patients with multiplicative perturbations on parameters $k_{10},k_{12}$ and $k_{21}$. Each parameter is multiplied by a constant drawn from a uniform distribution on the interval $\left[0.85,1.15\right]$. Moreover, the order of the fractional derivative $\alpha=1-\hat{p}/q$ is given by a random choice of $\hat{p}$ from a set of discrete values $\{17,18,20,21\}$, each one having the same probability. Denominator $q$ is fixed at $q=46$, while the nominal numerator is $p=19$. We simulate each patient by applying the same optimal drug scheduling sequence that was computed for the nominal one, that is, without online information about the state of each patient. Results are shown in Figure 10.

4 Conclusions

This paper gives an overview of the state of the art in numerical methods for the simulation of fractional-order systems along with validation results on a fractional-order pharmacokinetic model taken from the literature. Results are shown regarding the solution of an open-loop optimal control problem for the administration of Amiodarone to a patient whose pharmacokinetic parameters are assumed to be perfectly known. Additionally, we presented optimal control results for different patients whose pharmacokinetic parameters are not know perfectly. This way, we estimate and demonstrate how sensitive the administration scheduling is.