Eur J Drug Metab Pharmacokinet (2011) 36:189–196 DOI 10.1007/s13318-011-0057-6
ORIGINAL PAPER
A nonlinear two compartmental fractional derivative model Jovan K. Popovic´ • Diana Dolic´anin • Milan R. Rapaic´ • Stevan L. Popovic´ • Stevan Pilipovic´ • Teodor M. Atanackovic´
Received: 5 June 2011 / Accepted: 14 July 2011 / Published online: 30 July 2011 Ó Springer-Verlag France 2011
Abstract This study presents a new nonlinear two compartmental model and its application to the evaluation of valproic acid (VPA) pharmacokinetics in human volunteers after oral administration. We have used literature VPA concentrations. In the model, the integer order derivatives are replaced by derivatives of real order often called fractional order derivatives. Physically that means that the history (memory) of a biological process, realized as a transfer from one compartment to another, is taken into account with the mass balance conservation observed. Our contribution is the analysis of a specific nonlinear two compartmental model with the application in evaluation of J. K. Popovic´ (&) Department of Pharmacology, Toxicology and Clinical Pharmacology, Faculty of Medicine, University of Novi Sad, Hajduk Veljkova 3, 21000 Novi Sad, Serbia e-mail:
[email protected] D. Dolic´anin Department of Mathematics and Informatics, Faculty of Science and Mathematics, University of Novi Pazar, 36300 Novi Pazar, Serbia M. R. Rapaic´ Faculty of Technical Sciences, Computing and Control Department, University of Novi Sad, 21000 Novi Sad, Serbia S. L. Popovic´ Department of Hematology, Faculty of Medicine, University of Novi Sad, 21000 Novi Sad, Serbia S. Pilipovic´ Department of Mathematics and Informatics, Faculty of Science and Mathematics, University of Novi Sad, 21000 Novi Sad, Serbia T. M. Atanackovic´ Faculty of Technical Sciences, Institute of Mechanics, University of Novi Sad, 21000 Novi Sad, Serbia
VPA pharmacokinetics. The agreement of the values predicted by the proposed model with the values obtained through experiments is shown to be good. Thus, pharmacokinetics of VPA after oral application can be described well by a nonlinear two compartmental model with fractional derivatives of the same order proposed here. Parameters in the model are determined by the leastsquares method and the particle swarm optimization (PSO) numerical procedure is used. The results show that the nonlinear fractional order two compartmental model for VPA pharmacokinetics is superior in comparison to the classical (integer order) linear two compartmental model and to the linear fractional order two compartmental model. Keywords Fractional kinetics Nonlinear compartmental model Computer simulation Valproic acid
1 Introduction Fractional order models are found to be more adequate in many cases (Dokoumetzidis and Macheras 2009; Dokoumetzidis et al. 2010a; Kytarilos et al. 2010; Popovic et al. 2010a, b; Verotta 2010a, b), especially when a process under consideration exhibits hereditary behavior (Diethelm 2010). The interested reader can find detailed references regarding fractional calculus and fractional order models in Podlubny (1999), Hilfer (2000), Kilbas et al. (2006), and Diethelm (2010). Here, we simply review the basic concepts that will be utilized in the sequel. In our previous study, we have shown that a multidimensional fractional order compartmental model (Popovic et al. 2010a) gives the best agreement of the fitted curve to the experimental points if one uses the average value of
123
190
Eur J Drug Metab Pharmacokinet (2011) 36:189–196
a, denoted by a, significantly different from 1 although the case a ¼ 1 ð a 1Þ, can appear, as well. Since the individual approach to the patient is always advisable, it is completely understandable that for some patients a significantly differs from 1 ð a 6¼ 1Þ. Here we note that the appearance of integer order derivative models (with a = 1) is just a special case within the theory involving fractional order derivatives. So, in our approach, a is an individual parameter and can be used for individual optimization of the fitted curve to the individual concentrations with 0 \ a B 1. In our previous study (Popovic et al. 2010a), we determined a for slow-release diclofenac in healthy volunteers. In this study, we have two aims. First, we want to further improve and generalize mathematical model of fractional order. That is achieved by nonlinear approach, which is more general and contains the previous linear model, presented in our previous study (Popovic et al. 2010a), as a special case. Second, we want to check a for some other drugs. This is achieved with VPA concentrations in human volunteers for experiments presented in paper of Dutta et al. (2005).
2 Theoretical 2.1 Fractional calculus The integral of order a [ [0, 1] of a function f with respect to time variable t is defined as Zt 1 a f ðsÞðt sÞa1 ds; ð1Þ 0 It f ¼ CðaÞ 0
where C(a) is the Euler’s Gamma-function. Note that, for integer values of a = n, the above definition reduces to the well-known notion of the ordinary, n-fold integral. There are numerous ways to define the fractional derivative (Podlubny 1999; Hilfer 2000; Kilbas et al. 2006; Diethelm 2010). In the present paper, the Caputo definition is adopted (Diethelm 2010). Let us denote the ordinary n-fold derivative of a function f by f(n). The Caputo derivative of order a [ [n - 1, n] of a function f with respect to the time variable t is defined as Zt 1 a na ðnÞ f f ðnÞ ðsÞðt sÞna1 ds: ¼ 0 D t f ¼ 0 It Cðn aÞ 0
ð2Þ Although, it is not obvious from the above definition, one can show that for integer values of n, i.e., a = n we have 0 Dnt f ¼ f ðnÞ ðtÞ. Thus, the Caputo derivative of integer order reduces to the ordinary derivative.
123
The definition of the fractional integral and derivative is in a good accordance with the Laplace transform. The Laplace transform of an exponentially bounded function f is defined as Z1 f ðtÞest dt: ð3Þ ½Lf ðtÞðsÞ ¼ f^ðsÞ ¼ 0
The following facts hold for a [ (0, 1] ½L0 Dat f ðsÞ ¼ sa f^ðsÞ sa1 f ð0Þ;
ð4Þ
½L0 Ita f ðsÞ ¼ sa f^ðsÞ:
ð5Þ
It can also be proved that a a 0 It 0 Dt f
¼ f ðtÞ f ð0Þ:
ð6Þ
2.2 Nonlinear fractional two compartmental model Let us consider a two compartmental biological system after the paravascular drug application. The first compartment represents the place where the drug is applied, like the digestive tract, a muscle or the subcutaneous tissue. The second compartment represents the plasma and other body regions where the kinetic of the drug is uniform as in plasma. Conventionally, these two compartments are described by a system of ordinary differential equations (Popovic et al. 2010a). For example, q_ 1 ðtÞ ¼ k21 q1 ðtÞ;
ð7Þ
q_ 2 ðtÞ ¼ k21 q1 ðtÞ k02 q2 ðtÞ;
ð8Þ
where, q1 and q2 denote amounts of the drug in the first and second compartment, respectively. It has been proposed recently that fractional order models would better suit the dynamics of biological systems than the integer order ones. Among the first was the study of Dokoumetzidis and Macheras (2009). In our recent work (Popovic et al. 2010a), multi-compartmental models of fractional order were investigated. In the present study, we investigate a two compartmental model of the form sa11 1 0 Dat 1 q1 ðtÞ ¼ k21 f ðq1 ðtÞ; tÞ; sa22 1 0 Dat 2 q2 ðtÞ ¼ k21 f ðq1 ðtÞ; tÞ k02 q2 ðtÞ;
ð9Þ ð10Þ
where s1 and s2 are characteristic times for each compartment, a1 and a2 are constants satisfying 0 \ a1 B 1, 0 \ a2 B 1 and k21 and k02 are the constant factors specific to each individual. The initial amounts of the compartments are assumed to be: q1(0) = d1( = dose) and q2(0) = 0. In the analysis that follows we assume that a1 = a2 and s1 = s2 = s. Next, we discuss the conservation of matter problem for the proposed model. Namely, if there is no transfer rate between compartment 2 and
Eur J Drug Metab Pharmacokinet (2011) 36:189–196
external medium (i.e., k02 = 0) we must have the total amount of drug preserved. Thus, q1(t) ? q2(t) = q1(0) = const. By setting k02 = 0 and by adding Eqs. 9 and 10 we obtain 0 Dat 1 q1 ðtÞ þ 0 Dat 2 q2 ðtÞ ¼ 0 Now, since we assumed a1 = a2 = a, we have 0 Dat q1 ðtÞ þ0 Dat q2 ðtÞ ¼0 Dat ðq1 þ q2 Þ ¼ 0. Using the fact that Caputo derivative of a constant is equal to zero, we obtain q1(t) ? q2(t) = q1(0). Note that this relation is, in fact, a first integral of the system of differential Eqs 9 and 10. We note that the problem of conservation of matter was discussed recently in Dokoumetzidis et al. (2010b). For another approach to the conservation of mass problem in fractional calculus modeling (see Wheatcraft and Meerscharet 2008). As we stated in Popovic et al. (2010a), experimental results show that for all subjects the concentration in first compartment was not observed to rise immediately after the drug is administered, but after some time, i.e., there was a certain time lag after which the concentration starts to rise. Thus, for each subject first, a time instant tdelay was determined from the experimental data. It is defined by q1(t) = 0, q2(t) = 0, t [ [0, tdelay]. Then the system (9), (10) was solved in the time interval t [ [tdelay, T] with the initial conditions q1(tdelay) = d1, q2(tdelay) = 0. Thus, the final form of the equations that we will solve is a t 2 tdelay ; T ; ð11Þ 0 Dt q1 ðtÞ ¼ K21 f ðq1 ðtÞ; tÞ; a ð12Þ 0 Dt q2 ðtÞ ¼ K21 f ðq1 ðtÞ; tÞ K02 q2 ðtÞ; t 2 tdelay ; T ; with q1(tdelay) = d1, q2(tdelay) = 0. Here, we used K21 ¼ ; K02 ¼ ks02a : For specified f from the experimental data, we will determine the parameters a, K21, K02, tdelay and parameters involved in f. Several assumptions are adopted for the function f: k21 sa
1.A
2.A
3.A
4.A
f is a continuous function of both amount and time. This is because the drug amount is, by nature, a continuous quantity. f is non-negative at any time instance t and for any admissible drug amount q1. This is because f models the amount of drug transferred from the first to the second compartment. Negative f would mean that a certain amount of drug is traversing in the opposite direction. f is increasing function of the drug amount q1. This implies that higher drug concentration results in a higher transfer rate between the two compartments. It is also assumed that f(0, t) = 0 for any t. If this assumption would not be satisfied, then a certain amount of drug would be passed to the second compartment although the first one is empty. f is bounded and has bounded first derivative with respect to both arguments. These assumptions are mathematical prerequisites, used primarily to establish
191
existence and uniqueness of the solution, i.e., the well posedness of the proposed model. Further discussion is given in the Appendix. We used the function f in (11), (12) in the form b p q1 f ðq1 Þ ¼ 1 cos ; 2 d1
ð13Þ
where b is an arbitrary real number. We call this a nonlinear model. The form of f is chosen so that the speed of transfer (i.e., 0 Dat q1 ðtÞ) decreases with the amount of drug transferred from the first to the second compartment. Next, for comparison we used f in the form f ðq1 Þ ¼ q1 ;
ð14Þ
and call this case a linear model. Parameters for both models are determined through the least square numerical procedure. Remark One can formulate also a two compartmental model with nonlinearity in both equations, i.e., in Eq. 10 one may assume the nonlinear dependence of 0 Dat 2 q2 ðtÞ on q2. Also three compartment non-linear model similar to (11), (12) may be formulated. However, since in that case the analysis is more complicated and we do not use it in the experimental analysis, we do not present it in this paper.
3 Numerical procedure Solving fractional differential equations, especially nonlinear ones, is not an easy task. Therefore, only a brief outline of the solution procedure is presented here. A detailed account and analysis of the methods utilized in the solution process are given in the Appendix. The first nonlinear fractional equation in the proposed model (Eq. 9) was solved by transforming it to corresponding integral equation (see Eqs. 22, 23 in the Appendix). Then, resulting integral equation is solved using the successive approximations method, i.e. qkþ1 1 ðtÞ ¼ q0
k21 a k 0 It f ðt; q1 Þ; sa1 1
k 1;
ð15Þ
with q01 ðtÞ ¼ q0 ;
ð16Þ
where qk1 ðtÞ denotes the approximate solution obtained at the kth iteration. The actual solution is obtained in the limit, by allowing k tend to infinity. For numerical purposes, however, a finite number of iterations are used. In most cases, considered in the current study, the difference between the two consecutive approximations dropped below 0.001 after less than 500 iterations.
123
192
Eur J Drug Metab Pharmacokinet (2011) 36:189–196
Once the solution to the first equation, q1, is known, one can find and solve the second, linear equation (10) in the usual way, by means of the Laplace transform sa2 q^2 ðsÞ q2 ð0Þ ¼ k21 sa1 q^1 ðsÞ k02 q^2 ðsÞ:
ð17Þ
search procedure. In rare cases, additional 120 iterations were used to fine-tune the solution.
4 Experimental
Since q2(0) = 0, we obtain q^2 ðsÞ ¼
k21 sa1 q^1 ðsÞ : sa þ k02
ð18Þ
We used VPA concentrations in human volunteers published earlier (Dutta et al. 2005), where details about experimental setting are presented.
By applying the inverse Laplace transform it follows q2 ðtÞ ¼ ðL1 ðsa þ k02 Þ1 Þ ðk21 0 Ita f ðq1 ðtÞÞÞ;
ð19Þ 5 Results
where * denotes the convolution of functions. Solutions for q1(t) and q2(t) can be obtained from with appropriate initial data d1 and d2 (amounts of drug in compartments 1 and 2 at t = tdelay). We analyze the concentration c2 = q2(t)/v2 in the second compartment. Since concentrations depend on the volumes of compartments, we consider quotient d1/v2 = d as a parameter which depends on each subject separately. In the current study, all measurements were scaled so that concentration values are expressed in multiples of 100 mg/cm3. The time is also normalized, expressed in multiples of 8 h. The fractional integral was calculated using a fixed step length discrete approximation of the integration process. The same is true for the convolution operator. The time grid used in numerical simulations was created using D = 0.01 as the sample time (corresponding to 4.8 min in physical time). Delay is implemented as an appropriate time-shift of the calculated concentrations. Further details are given in the Appendix. For each of the experimental data set, five parameters were estimated for linear models (Popovic et al. 2010a) and 1 more for the nonlinear ones: k21, k02, a, d and tdelay, and also b for the nonlinear model. For a particular choice of these parameters, the mean-square error (MSE) of the concentration estimate c2(t) is MSE ¼
M 1X ½c2 ðti Þ cm ðti Þ2 ; M i¼0
ð20Þ
where cm(ti) is the measured concentration, ti are times of measurements and M is the number of measurements. It is clear that MSE is a function of the model parameters. The optimal model parameters are those for which the MSE value is minimal. Since parameters enter in a highly nonlinear, non-convex way, particle swarm optimization (PSO) algorithm is used during the optimization process (Rapaic´ 2010). Like some other biologically inspired optimizers, PSO is known for its robustness and the ability to find good solutions of difficult optimization problems. Most of the examples in the current study are obtained using PSO with single swarm of 25 particles and 120 iterations of the
123
In the sequel, we analyze experimental data to confirm the flexibilities of our general concept. The results are shown in Fig. 1 for the mean plasma VPA concentration–time profile for fast release formulation in Dutta et al. (2005). The values of parameters and MSE for linear fractional order model (11), (12), (14) and for nonlinear fractional order model (11), (12), (13) are presented in Table 1. As could be seen, the agreement of experimental values with the values obtained for nonlinear model is good. In Fig. 2, we show results for the mean plasma VPA concentration–time profile for medium release formulation in Dutta et al. (2005). The values of parameters and MSE for linear fractional order model (11), (12), (14) and for nonlinear fractional order model (11), (12), (13) are presented in Table 1. Now the nonlinear model shows better agreement with experimental data than linear. In Fig. 3, we show results for the mean plasma VPA concentration–time profile for slow release formulation in Dutta et al. (2005). The values of parameters and MSE for
Fig. 1 Measured VPA concentrations (symbols) reported in literature (Dutta et al. 2005) for fast release formulation and the predicated (linear and nonlinear fractional order models) values (lines)
Eur J Drug Metab Pharmacokinet (2011) 36:189–196
193
Table 1 Parameter estimates after fitting Eq. 12 (linear model) and Eq. 13 (nonlinear model) to the literature data (Dutta et al. 2005) Formulation Fast release Medium release Slow release
Model
a
k21
k02
d
tdelay
b
MSE 0.000103120
Linear
0.887399
3.58349
2.12001
0.862407
0.0589739
–
Nonlinear
0.832752
2.65534
5.26016
2.312580
0.0622120
3.87584
0.000159142
Linear
1
1.88990
1.26536
0.570275
0.0189203
–
0.000195894
Nonlinear
0.990236
0.93174
2.92929
1.289710
0.0090271
6.91051
0.000027852
Linear
1
1.65762
1.07898
0.479521
0.0014876
–
0.000538844
Nonlinear
0.894522
0.61617
2.08185
0.769679
0.0132396
139.909
0.000162537
linear model (11), (12), (14) and for nonlinear model (11), (12), (13) are shown in Table 1. Again the nonlinear model shows better agreement with experimental results.
6 Discussion Our previous applications of linear fractional order two compartmental model on diclofenac pharmacokinetics (Popovic et al. 2010a, b) are an extension of the results presented by Dokoumetzidis and Macheras (2009), with fractional order one-compartmental model. Our general approach to linear fractional order two compartmental models (Popovic et al. 2010a) contains the classical case integer order model (a = a1 = a2 = 1) as well the case a = a1 = a2 [ 1, and a = a1 = a2 \ 1. In this study, as in the paper of Kytarilos et al. (2010), experimental VPA data of Dutta et al. (2005) were used to examine the method. But in Kytarilos et al. (2010) the problem was treated mathematically as a one-compartment model.
Fig. 3 Measured VPA concentrations (symbols) reported in literature (Dutta et al. 2005) for slow release formulation and the predicated (linear and nonlinear fractional order models) values (lines)
By a nonlinear approach which contains the previous linear fractional order model, presented in Popovic et al. (2010a) as a special case, for f(q1) = q1, we further generalize and improve mathematical model of fractional order as seen in Figs. 2 and 3. From our results, we learned that kinetics of VPA is of a slow type (Popovic et al. 2010a) (when a \ 1) or standard one (when a = 1). Note that d is a parameter depending on v2, and because of that it is different for every subject.
7 Conclusions Two major conclusions may be drawn from our analysis. The first one is that the nonlinear fractional compartmental method could be successfully used to model pharmacokinetics. The second one is that the kinetics of VPA corresponds to the order of derivative a less than or equal to one. Fig. 2 Measured VPA concentrations (symbols) reported in literature (Dutta et al. 2005) for medium release formulation and the predicated (linear and nonlinear fractional order models) values (lines)
Acknowledgments This research is supported by Republic of Serbia, Ministry of Science grants 174024 (SP), 174005 (TA) and 32018 (ZJ) and Autonomous Province of Vojvodina, Provoncial Secretariat
123
194
Eur J Drug Metab Pharmacokinet (2011) 36:189–196
for science and technological development grant 114-451-2048/201101 (JP).
Appendix In the sequel, we present a formal result regarding existence and uniqueness for the equations of the proposed model. The results are basis for the iterative procedure that can be used to construct the solution to a nonlinear fractional initial value problem (11), (12). We also propose a numerical procedure for evaluation of fractional order integrals and give an error estimate for such a procedure. The procedure is finally used to reduce a linear fractional differential equation to a system of linear algebraic equations, which can be solved in a standard way.
One can always choose t1 [ (0, T] so that t1a L\1; a CðaÞ
and thus, according to the Banach fixed point theorem, the unique solution to the initial value problem is qðtÞ ¼ limk!1 qk ðtÞ: The solution constructed in such a way is defined only for t [ [0, t1]. However, this solution can be extended on the interval t [ [t1, 2t1]. For such a t, (23) becomes qðtÞ ¼ q0 þ 0 Ita Fðt; qðtÞÞ Zt1 1 ¼ q0 þ Fðt; qðtÞÞðt sÞa1 ds Cða 1Þ 0
¼ q0 þ Reduction to integral equation
Theorem 1. Let a [ (0, 1) and let F:[0, T] 9 D, D (
