- Email: [email protected]
Two-level time-domain decomposition based distributed method for numerical solutions of pharmacokinetic models
Computers in Biology and Medicine 41 (2011) 221–227
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Two-level time-domain decomposition based distributed method for numerical solutions of pharmacokinetic models Li Liu a,c,n, Choi-Hong Lai b, Shao-Dan Zhou c, Fen Xie c, Rui Lua a
School of IoT Engineering, Jiangnan University, Wuxi 214122, Jiangsu Province, China School of Computing and Mathematical Sciences, University of Greenwich, Old Royal Naval College, Greenwich, London SE10 9LS, UK c Suzhou University Affiliated Fourth People’s Hospital, Wuxi Fourth People’s Hospital, Wuxi 214062, Jiangsu Province, China b
a r t i c l e i n f o
abstract
Article history: Received 29 March 2010 Accepted 4 February 2011
In order to predict variations of drug concentration during a given period of time, numerical solutions of pharmacokinetic models need to be obtained efficiently. Analytical solutions of linear pharmacokinetic models are usually obtained using the Laplace transform and inverse Laplace tables. The derivations of solutions to complex nonlinear models are tedious, and such solution process may be difficult to implement as a robust software code. For nonlinear models, the fourth-order Runge–Kutta (RK4) is the most classical numerical method in obtaining approximate numerical solutions, which is impossible to be implemented in distributed computing environments without much modification. The reason is that numerical solutions obtained by using RK4 can only be computed in sequential time steps. In this paper, time-domain decomposition methods are adapted for nonlinear pharmacokinetic models. The numerical Inverse Laplace method for PharmacoKinetic models (ILPK) is implemented to solve pharmacokinetic models with iterative inverse Laplace transform in each time interval. The distributed ILPK algorithm, which is based on a two-level time-domain decomposition concept, is proposed to improve its efficiency. Solutions on the coarser temporal mesh at the top level are obtained one by one, and then those on the finer temporal mesh at the bottom level are calculated concurrently by using those initial solutions that have been obtained at the top level decomposition. Accuracy and efficiency of the proposed algorithm and its distributed equivalent are investigated by using several test models. Results indicate that the ILPK algorithm and its distributed equivalent are good candidates for both linear and nonlinear pharmacokinetic models. Crown Copyright & 2011 Published by Elsevier Ltd. All rights reserved.
Keywords: Pharmacokinetics Distributed computing Numerical inverse Laplace Runge-Kutta
1. Introduction Pharmacokinetic models are used to describe the absorption, distribution, metabolization and elimination (ADME) of drugs in vitro, in loci and in vivo [1]. In order to predict variations of drug concentrations during a given period of time, numerical solutions of the models need to be obtained efficiently. For linear models, the Laplace transform method supplemented with a lookup inverse Laplace table is often used to obtain the corresponding analytical solutions [2]. For nonlinear models, the fourth-order Runge–Kutta (RK4) [3] is a typical numerical method widely used to obtain approximate solutions, while the Adomian decomposition method has been reported to be used to obtain semianalytical solutions [4]. In practice, derivations of analytical and semi-analytical solutions are tedious and such solution process is
n Corresponding author at: School of IoT Engineering, Jiangnan University, Wuxi 214122, Jiangsu Province, China. E-mail address: [email protected] (L. Liu).
also difficult to be implemented as a robust software code. On the other hand, numerical solutions along the temporal axis calculated by using the RK4 method can only be obtained in sequential time steps. Inevitably, computing time of such temporal integration method becomes significant for nonlinear problems. The goal of this paper is to present an alternative approach to obtain numerical solutions of both linear and nonlinear pharmacokinetic models. At the same time, the efficiency of the method implemented using the distribution strategy is discussed and tested. Our approach is based on a time-domain decomposition method. There are a number of time-domain decomposition methods (TDDM). The key point of the time-domain decomposition is to split up the temporal axis into smaller intervals, which are used for the transmission of data across the intervals [5]. In this paper, a time-domain decomposition method is applied to solve a first-order nonlinear system. The ILPK algorithm is implemented to deal with the iterative inverse Laplace method in each time interval. One property that domain decomposition methods have is the suitability for parallel computing [6]. In our previous work, a two-level time-domain decomposition method
0010-4825/$ - see front matter Crown Copyright & 2011 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2011.02.003
222
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
was applied to obtain numerical solutions of time-dependent nonlinear problems for European options [7]. In this paper, it is used to decouple the numerical calculation process of pharmacokinetic models. At the top level, the entire time-domain is divided into coarser time intervals. At the bottom level, each coarser time interval is further sub-divided into finer time intervals. Solutions on the coarser temporal mesh are obtained one by one in a sequential manner, and those on the finer temporal mesh are calculated concurrently with initial solutions that have been obtained at the top level decomposition. Further distributed properties are also induced into the algorithm at the coarser temporal computation by using Laplace transformation. In the numerical experiments, the sequential and the distributed methods are implemented and tested for several models chosen from pharmacokinetics and both the accuracy and the efficiency of the proposed methods are investigated. Parameter settings of these methods are also discussed.
problem of solving linear or nonlinear model is converted into obtaining solutions of Eq. (3). In the next section, numerical methods are given in details.
3. Numerical methods to solve pharmacokinetic models 3.1. The Runge–Kutta methods The Runge–Kutta methods represent a class of numerical methods in obtaining approximate numerical solutions of differential equations related to initial value problems [11]. The methods rely on the Taylor series expansion but do not require the computation of high order derivatives of the right hand side function of the differential equations. The fourth-order Runge– Kutta (RK4) is a typical numerical method widely used to obtain approximate solutions, which requires four evaluations of updates per time step. Its error per time step Dt is of the order of Dt5, while the total accumulated error is of the order of Dt4.
2. The pharmacokinetic models Compartment model [8] is a classical type of pharmacokinetic models and views the body as a system composed of many compartments taking into variations of drug concentration. The distribution of a drug is described usually as transportation of drug concentration between compartments and is represented mathematically by a set of differential equations. For example, an n-compartment model in Ref. [9] is defined as follows: x_ i ðtÞ ¼
n X
kij xi ðtÞ þ
j¼1 iaj
n X
kji xj ðtÞ,
i,j ¼ 1,2,3,. . .,n
ð1Þ
j¼1 i aj
with initial conditions xi(0)Z0. Here xi(t) represents the drug concentration in the ith compartment at time t, kij is the proportional constant associated with the rate of absorption or elimination from the ith compartment to the jth compartment. In addition, there are many drugs which are eliminated by enzymes and transporters expressed mainly in the intestine, liver and kidney, and are bound to proteins in the blood and tissues. Since the numbers of enzymes, transporters and binding proteins are limited, the disposition of drugs is capacity-limited. In this context, the Michaelis–Menten model is often used to describe the rapid equilibrium kinetics of enzyme systems [10]. It is also embedded into compartment models to represent the rate of change from one compartment to another. Here is an example of one compartment model with Michaelis–Menten elimination. The drug is given by intravenous injection. Let Vm be the maximum elimination rate, and Km be the Michaelis–Menten constant. The coefficient Vm/(Km + x(t)) represents the capacity-limited clearance. The variation of drug concentration described by Michaelis–Menten equation is defined as below _ ¼ xðtÞ
Vm xðtÞ Km þ xðtÞ
ð2Þ
A generalized mathematical description of pharmacokinetic models capable of describing both linear and nonlinear compartment models can be written as X_ ¼ AX þHðXÞ þ FðtÞ, t A ½0,T
ð3Þ
where X is the concentration matrix of compartments and X_ ¼ dX=dt. A multi-compartment model composes three parts: the linear part AX represents first-order kinetic distribution between compartments where A is a coefficient matrix, the nonlinear part H(X) representing nonlinear absorption or elimination and the external input function F(t) which is used to express the administration protocol of drugs. As a result, the
3.2. Numerical inverse Laplace and time-domain decomposition The Laplace transform of the time dependent function x(t) is an integral transform defined by the integral Z 1 elt xðtÞdt ¼ UðlÞ ð4Þ LðxðtÞÞ 0
It is particularly useful in solving linear ordinary differential equations where an ordinary differential equation with temporal derivatives is transformed into a form that does not involve time variable. In this paper the authors examine an application of Laplace transformation to nonlinear differential equations. There are many different methods for numerically inverting the Laplace transform. Among them, the Graver–Stehfest inverse Laplace method which was proposed by Graver [12] and derived by Stehfest [13] gave good accuracy on a fairly wide range of problems [14]. It should be noted that the Laplace integral is known to exist on the condition that xðtÞ be piecewise continuous on every finite interval in t Z 0 and satisfy xðtÞ r Meat for certain constants M and a. It is clear that drug concentrations are under this condition. Hence, taking Laplace transform of Eq. (3) leads to ð5Þ UðlÞ ¼ ½lIA1 ½LðHðXÞÞ þLðFðtÞÞ þ Xðt0 Þ subject to the constraint lIA a0. Here l A flj ,j ¼ 1,2,. . .,mg is a finite set of transformation parameters defined as below
lj ¼ j
ln 2 , T
j ¼ 1,2,. . .,m
ð6Þ
In order to retrieve X, the numerical inverse Laplace transform due to Graver–Stehfest method as listed in Eq. (7) may be used X
m ln 2 X U w UUðlj Þ T j¼1 j
ð7Þ
where the weighting factor wj is given by 2 3 minðj,m=2Þ X im=2 ð2iÞ! j þ ðm=2Þ 4 5 U wj ¼ ð1Þ ðm=2iÞ!i!ði1Þ!ðjiÞ!ð2ijÞ! i ¼ ð1 þ jÞ=2 In the Graver–Stehfest method, m is the only parameter to be set. As defined in Eq. (7), the approximation of Laplace transformed X is obtained by integrating m different parametric functions of Uðlj Þ. The higher the value of m, the more computing time is required to calculate X but not necessarily more accurate. Note that it is possible to compute these parametric functions simultaneously using a distributed computing environment.
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
The reason is that these parametric functions are independent from each other. To solve the nonlinear part of the model, time-domain decomposition is applied firstly. The key point of the time-domain decomposition is to split up the temporal axis ½0,T into smaller intervals Ip ¼ ½tp ,tp þ 1 with t0 ¼ 0, tN ¼ T, tp otp þ 1 , and p ¼ 0,1,. . .,N1, in order to transmit data across the boundary of Ip . As such the numerical inverse Laplace method has to be done sequentially in each time interval. As in Eq. (8), HðXÞ is computed using the linearization X~ , which is updated within an iterative process until convergence is achieved for HðXÞ within that time interval. Each step of the iterative update process involves a numerical solution to the equation X_ ðtÞ ¼ AX þ HðX~ Þ þ FðtÞ,
t A ½tp ,tp þ 1
ð8Þ
Let Xðtp Þ and Xðtp þ 1 Þ be the numerical solutions of Eq. (3) at tp and tp þ 1 . A Laplace transform is applied to Eq. (8), being defined on t A ½tp ,tp þ 1 , in its differential form which leads to UðX~ , lÞ ¼ ½lIA1 ½HðX~ Þ=l þ LðFðtÞÞ þXðtp Þ. In this paper, the iterative inverse Laplace to obtain numerical solution Xðtp þ 1 Þ of nonlinear pharmacokinetic models, using Xðtp Þ as the initial approximation to X~ is implemented. The approximation X~ of Xðtp þ 1 Þ is updated with X ðkÞ ðtp þ 1 Þ at the kth step iteration until the convergence criteria :X ðkÞ ðtp þ 1 ÞX ðk1Þ ðtp þ 1 Þ: o e is achieved. Here is the main function the ILPK algorithm used to compute the solution X. FUNCTION X(tp + 1)¼ILPK(X(tp)) Initial approximation: X(0)(tp + 1) ¼X(tp); k¼0; // Iterative Inverse Laplace in Ip ¼[tp,tp + 1] ITERATE X~ ¼ X ðkÞ ðtp þ 1 Þ; Compute simultaneous for j¼ 1, 2, y, m: UðX~ , lj Þ ¼ ½lj IA1 ½HðX~ Þ=lj þ LðFðtÞÞ þ Xðtp Þ, where lj ¼ j tp þln1 2tp ; Compute X ðk þ 1Þ ðtp þ 1 Þ
ln 2 tp þ 1 tp
" j þ ðm=2Þ
where wj ¼ ð1Þ
U
m P
wj UðX~ , lj Þ,
j¼1
minðj,m=2Þ P i ¼ ð1 þ jÞ=2
im=2 ð2iÞ! ðm=2iÞ!i!ði1Þ!ðjiÞ!ð2ijÞ!
# ;
k ¼k+1; UNTIL :X ðkÞ ðtp þ 1 ÞX ðk1Þ ðtp þ 1 Þ:o e; RETURN X(k)(tp + 1); END-FUNCTION
If H(X) ¼0, Eq. (3) represents a linear model, and its solution at time t ¼tp + 1, p ¼0, 1, ..., N 1 can be obtained by X(tp + 1)¼ ILPK(X(t0)) simultaneously when distributed to N processors. However, for the nonlinear model, solutions can only be obtained sequentially since computing X(tp + 1) is very depended on its initial solution X(tp). Section 4 proposes a method of improving the efficiency of the ILPK algorithm under distributed environments—two-level time-domain decomposition.
4. Distributed algorithms based on two-level time-domain decomposition In our previous work, a two-level time-domain decomposition method was used to decouple time-dependent nonlinear problems for European options. The main idea of the two-level timedomain decomposition is to divide the entire time-domain [0,T] into smaller time intervals. At the bottom level, [0,T] is equally divided into finer time intervals Ibp ¼[tp,tp + 1], with, tp + 1 ¼tp + Dt, p ¼0, 1, ..., N 1, t0 ¼0, tN ¼T. At the top level, Nfine pieces of finer
223
Fig. 1. Illustration of the process of the distributed ILPK based on the two-level time-domain decomposition with DT¼ 5Dt
time intervals are integrated into a coarse time interval, Itq ¼ [Tq,Tq + 1], with Tq + 1 ¼Tq +Nfine Dt, q¼0, 1, ..., M 1, M ¼ N=Nfine and T0 ¼t0. Based on the two-level time-domain decomposition, X(Tq), q¼0, 1, ..., M 1 on the coarser temporal mesh are computed one by one in a sequential manner, and are being used as initial conditions of solutions on the finer time intervals. So remaining solutions on the finer intervals X(Tq + ti),ti ¼i Dt, and those X(Tk + ti), k aq, i¼ 1, ..., Nfine 1 can be obtained concurrently. In short, sequential X(tp) is mapped into two-level time domain X(Tq + ti), q ¼ p=Nfine , i ¼ kmodNfine . Given DT¼5Dt, the process of the distributed ILPK based on the two-level timedomain decomposition is illustrated in Fig. 1. For example, after X(t5), X(t10), X(t15) have been calculated sequentially, X(t1) is obtained with initial condition X(t0) followed by X(t2), X(t3), X(t4) and parallelized with X(t6), X(t11)y The pseudo-code of the distributed ILPK algorithm is described as below. FUNCTION Distributed ILPK IF H(X)¼0, // for linear models, do one-level time-domain decomposition THEN Distribute p ¼1, y, N tp to processors, X(tp) ¼ILPK(X(t0)); END-Distribute and receive X(tp). ELSE // for nonlinear models, do two-level time-domain decomposition T0 ¼t0; DT¼Nfine Dt; M ¼ N=Nfine ; FOR q¼0 To M 1, do coarser temporal mesh on tA[Tq,Tq + 1] sequentially Tq + 1 ¼Tq + DT; // top level time interval X(Tq + 1) ¼ILPK(X(Tq)); END-FOR Distribute q¼0 To M 1,X(Tq) to processors, For i¼1 To Nfine 1, do finer temporal mesh on tA[Tq + ti 1,Tq + ti] ti ¼i Dt; // bottom level time interval X(Tq + ti)¼ ILPK(X(Tq + ti 1)); END-FOR END-Distribute and receive X(Tq + t1), y, X(Tq + tNfine 1 ). END-IF END-FUNCTION
From the pseudo-code above we can see that if we replace the ILPK algorithm with the RK4, the computing of solutions can also be distributed based on the two-level time-domain decomposition method. In the next section, numerical experiments are provided to examine the accuracy and the efficiency of the distributed ILPK and RK4 algorithms.
5. Experiments and discussion In our experiments, the ILPK and the RK4 methods are programmed with MatlabTM software. In Section 5.1, linear,
224
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
nonlinear pharmacokinetic models, as well as an integrated pharmacokinetic–pharmacodynamic model are solved by using the ILPK algorithm. In Section 5.2, the accuracy of the distributed algorithms with variation of two-level time-domain distributed strategy DT ¼Nfine Dt is investigated. Theoretical and experimental efficiencies are calculated and observed with Matlab distributed computing platform in Section 5.3. Finally, parameter setting of the distributed algorithms is discussed in Section 5.4.
Solutions X(ILPK) and X(RK4) in tA[0,0.05] and their discrepancies are shown in Table 2. The discrepancies :XðILPKÞXðRK4Þ: are less than 4E-6.
5.1. Accuracy analysis of the ILPK algorithm
A more complicated model was a pharmacokinetic– pharmacodynamic (PK–PD) model, which was built by Barbolosi [17] to predict cancer chemotherapy. Here, we extend the application of the ILPK method to solve this model. c1(t) is the concentration of administered anticancer drug in the plasma. The mass balance for the two-compartment model leads to 8 1 ðtÞ < dcdt ¼ ðkE1 þ k12 ÞUc1 ðtÞ þ uðtÞ V1 ð11Þ : dc2 ðtÞ ¼ k12 U V1 Uc1 ðtÞkE2 Uc2 ðtÞ V2 dt
Example 1. A linear model for extra vascular administration A linear pharmacokinetic model for oral administration proposed in literature [15] is rewritten as Eq. (3) and displayed in Eq. (9). Here, H(X)¼0, F(t)¼0. Given initial dose, drug went into " # x1 compartments with first-order absorption process. X ¼ x2 represented the concentrations of the two compartments. The rate constants of the differential equation were k12 ¼2.008 and k21 ¼0.1855. The initial drug concentrations in the two compartments at t ¼0 were x1(0) ¼103.8613 and x2(0) ¼0. " # k12 0 _ X¼ X ð9Þ k21 k12 The classical method to obtain analytical solutions x1, x2 of the linear compartment model is the Laplace transform supplemented with a look-up inverse Laplace table. Here, analytical solutions of x1, x2 were deduced. The parameter of both the ILPK and RK4 methods were set as m¼ 16, and time step Dt ¼0.001. Approximate solutions obtained by the ILPK and RK4 method were denoted as x1(ILPK), x2(ILPK), x1(RK4) and x2(RK4). Discrepancies between analytical solutions and those obtained by the ILPK and the RK4 method :x1 x1 ðRK4Þ:, :x2 x2 ðRK4Þ:, :x1 x1 ðILPKÞ: and :x2 x2 ðILPKÞ: are recorded in Table 1. It is clear that the approximate solutions obtained by the ILPK algorithm are closer to analytical solutions than the RK4 method. Example 2. Nonlinear elimination
pharmacokinetic
model
for
drug
A model obeying the first order and Michaelis–Menten elimination kinetics in the case of extra vascular administration was displayed in Eq. (10), which was solved by the RK4 method [16]. Here, the model with constants ka ¼0.5 h 1, Vd/F¼40 L, 1 D ¼150 mg, tlag ¼0, Vm ¼ 0:8 mg L1 h , Km ¼10 mg L 1 is solved by using the ILPK method. In the numerical experiments the parameters of the ILPK algorithm were set as m¼10, Dt ¼0.0001.
dX Vm X FD ka ðttlag Þ ¼ kX þ ka e dt Km þX Vd
ð10Þ
Example 3. Pharmacokinetic–pharmacodynamic model of cancer chemotherapy
where 8 > < 0:205 uðtÞ ¼ 0:1 > : 0
for
0 rt o 1
for
1r t o22
for
22 rt
is the drug delivery by intravenous infusion. Assuming that the growth of tumor behaved according to a Gompertz-type growth equation, n(t) denotes the number of tumor cells at time t, and H(.) is a Heaviside step function dnðtÞ ¼ lUnðtÞUln y=nðtÞ kU½c2 ðtÞCMIN UnðtÞUHðc2 ðtÞCMIN Þ dt
ð12Þ
Table 2 Solutions of Example 2 obtained using the ILPK and RK4 methods. t(h)
X(ILPK)
X(RK4)
:XðILPKÞXðRK4Þ:
0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050
0.0093613 0.0186954 0.0280025 0.0372826 0.0465358 0.0557622 0.0649618 0.0741347 0.0832810 0.0924008
0.0093610 0.0186960 0.0280030 0.0372830 0.0465370 0.0557630 0.0649630 0.0741360 0.0832830 0.0924030
3.53E-7 7.05E-7 1.06E-6 1.41E-6 1.75E-6 2.10E-6 2.45E-6 2.79E-6 3.14E-6 3.48E-6
Notes: X(ILPK), X(RK4) were approximate solutions obtained using the ILPK and RK4 methods, respectively. :XðILPKÞXðRK4Þ: were discrepancies between these solutions.
Table 1 Discrepancies between analytical solutions and those obtained by the ILPK and RK4 methods. t(h)
x1
x2
:x1 x1 ðRK4Þ:
:x2 x2 ðRK4Þ:
:x1 x1 ðILPKÞ:
:x2 x2 ðILPKÞ:
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
84.7960 69.3696 56.7496 46.4255 37.9796 31.0702 25.4178 20.7937 17.0108 13.9161
18.8818 33.8143 45.6926 55.0786 62.4318 68.1279 72.4744 75.7224 78.0775 79.7078
0.1704 0.1394 0.1141 0.0933 0.0763 0.0625 0.0511 0.0418 0.0342 0.0280
0.1670 0.1332 0.1056 0.0831 0.0648 0.0498 0.0376 0.0278 0.0197 0.0132
3E-06 2E-06 1E-05 4E-06 1.1E-05 1.9E-05 2E-06 4.2E-05 8.8E-05 1.7E-04
0 1E-06 1E-06 4E-06 9E-06 1.2E-05 0 3.8E-05 1.1E-04 1.8E-04
Notes: x1, x2 were analytical solutions of Example 1. :x1 x1 ðRK4Þ:, :x2 x2 ðRK4Þ:, :x1 x1 ðILPKÞ:, :x2 x2 ðILPKÞ: were discrepancies between analytical solutions and those obtained by the ILPK and RK4 methods.
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
225
In addition, dynamical toxicities induced by drug exposures is also taken into account in a cumulative process expressing myelosuppression in the bone marrow site, which is formulated as dwðtÞ ¼ rc vUwðtÞmUwðtÞUc1 ðttÞ dt
ð13Þ
In the numerical experiments the constants of the model were V1 ¼ 20 L, kE1 ¼ 0.48 day 1, k12 ¼0.02 day 1, V2 ¼3 L, kE2 ¼0.2 day 1, l ¼3 10 3 day 1, y ¼1012, k ¼ 10 g1 L day1 , CMIN ¼1 mg mL 1, rc ¼0.8 10 3 L 1 day 1, v¼0.1 day 1, m ¼50 g 1 L day 1, t ¼10 days. The initial values of the model were c1(0)¼0 g L 1, c2(0)¼ 0 g L 1, n(0)¼n0, w(0)¼w0, no ¼5 1010 cells, wo ¼rc/v, tA [0,45] days. To present the variations of drug concentrations, the number of tumor cells, and the number of white blood cells during 45 days, the ILPK and the RK4 method were used. Since there was time delay for toxicities, Eq. (13) was solved after Eqs. (11) and (12) had been done. The parameters of the ILPK algorithm were set as m¼10, and Dt ¼0.005. Discrepancies between approximate solutions obtained by the ILPK and RK4 methods are recorded in Table 3. It indicates that solutions obtained by the ILPK algorithm are very closer to those of the RK4 method.
Fig. 2. Comparison of discrepancies between solutions of Example 3 obtained using the distributed ILPK, RK4 algorithms and using their sequential ones.
those of the distributed ILPK algorithm. It also indicates that the accuracy of the distributed ILPK algorithm is little affected when DT¼5Dt and 10Dt. 5.3. Efficiency analysis of the distributed algorithms
5.2. Accuracy analysis of the distributed algorithms It is clear that the time interval Dt¼tp + 1 tp is one main contributing factor to the accuracy of solutions obtained by both the sequential ILPK and the RK4 methods, since solution X(tp + 1) is obtained by using X(tp) as the initial approximation. Compared to its sequential equivalent, the distributed ILPK and RK4 methods calculate solutions in two phases. Firstly, solutions at the bottom level (finer intervals) are dependent on that at the top level (coarse intervals) solutions. It can be imagined that with a larger coarse interval DT¼Nfine Dt, discrepancies between solutions of the distributed algorithm and its sequential equivalent become greater. In this paper we investigate how far the accuracy of the distributed algorithms is being affected by the distributed strategy represented as coarse interval DT¼ Nfine Dt. The distributed ILPK and the distributed RK4 algorithms were run to solve Example 3 with DT¼ 5Dt,10Dt,25Dt,50Dt, respectively. The parameters of the algorithms used in the numerical experiments were m¼10 and Dt¼0.01. Accumulated discrepancies for 45 days between solutions obtained by using the distributed and its sequential algorithm are illustrated in Fig. 2. When DT increases from 5Dt to 50Dt, discrepancies between solutions obtained by the distributed algorithms and their sequential ones goes up quickly. In addition, the discrepancies between solutions of the distributed RK4 and its sequential version increased a little more quickly than Table 3 Solutions of Example 3 obtained by using the ILPK and RK4 methods. t(day)
c1(ILPK)
c2(ILPK)
n(ILPK)
w(ILPK)
d(c1)
d(c2)
d(n)
d(w)
5 10 15 20 25 30 35 40 45
0.00970 0.00998 0.01000 0.01000 0.00218 0.00018 0.00001 0.00000 0.00000
0.00361 0.00550 0.00624 0.00651 0.00501 0.00212 0.00080 0.00030 0.00011
9.52 7.92 6.22 4.80 3.75 3.38 3.37 3.43 3.48
0.08000 0.08000 0.01945 0.01370 0.01336 0.01333 0.02094 0.04034 0.05545
7E-05 2E-06 4E-06 3E-06 1E-04 1E-05 1E-05 0 0
7E-05 4E-05 1E-05 3E-06 8E-05 4E-05 2E-05 0 0
0.02 0.03 0.03 0.03 0.03 0.04 0.04 0.05 0.04
2E-06 1E-06 3E-04 2E-05 1E-05 1E-06 4E-04 4E-04 3E-04
Notes: c1(ILPK), c2(ILPK), n(ILPK), w(ILPK) were approximate solutions of Example 3 obtained by the ILPK method. dðc1 Þ ¼ :c1 ðILPKÞc1 ðRK4Þ:, dðc2 Þ ¼ :c2 ðILPKÞ c2 ðRK4Þ:, dðnÞ ¼ :nðILPKÞnðRK4Þ:, dðwÞ ¼ :wðILPKÞwðRK4Þ: were discrepancies between approximate solutions of Example 3 obtained using the ILPK and RK4 methods.
The execution time of a distributed algorithm composes of computing time and communication overheads. The latter is governed by the computing network environments. When problems to be handled are more complex and can be parallelized without much communication between different problems, the communication overheads can be neglected. Defining one work unit as the computational work required for solving the pharmacokinetic model X(tp + 1)¼ ILPK(X(tp)) in a given time interval [tp,tp + 1]. The total parallel work unit of the distributed ILPK algorithm is composed of the total work unit at the top level and the distributed work unit at the bottom level executed in each processor. In this paper, the theoretical efficiency of the distributed algorithm is defined as Eq. (14), with the number of parallel work unit NWU(parallel) divided by the number of sequential work unit NWU(sequential) it required. The smaller values of efficiency show that the algorithm is more efficient. It gives us a measure to evaluate how much time saving is achieved when using the distributed strategy without the effect of communication overheads. Theoretical efficiency ¼
NWUðparallelÞ NWUðsequentialÞ
ð14Þ
For linear models, X(tp + 1)¼ILPK(X(tp)), p¼0, 1, ..., N 1, the efficiency can be calculated simultaneously in processors. Given the distributed algorithm is executed in Np processors, its theoretical efficiency to solve linear model is 1/Np. When solving nonlinear problems with two-level time-domain decomposition DT¼Nfine Dt, the NWU(parallel) is composed of n=Nfine work units of initial conditions and those on finer time intervals in each processor n (Nfine 1)/(NfineNp). So, the theoretical efficiency of the distributed ILPK algorithm for solving nonlinear problem is efficiency¼1/Nfine + 1/Np 1/(NfineNp), which is inversely proportional to the number of processors Np and the coefficient of coarser time interval Nfine. Table 4 shows the theoretical efficiency of the distributed ILPK calculated with a range of coarse time intervals DT¼5Dt, 10Dt,25Dt,50Dt, and with different number of processors Np ¼2,4,8,10,16 for the problem described in Example 2. We can see that the distributed algorithm saves much more computing time with regard to more processors and bigger time intervals. In this paper, experimental efficiency with various Nfine, Np was also observed. The distributed environment was set up with MatlabTM Distributed Computing Engine (MDCE). The performance analysis was obtained on a PC cluster consisting of 16
226
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
Table 4 Theoretical efficiencies of the distributed ILPK algorithm to solve Example 2. Np
DT¼5Dt (%)
DT¼ 10Dt (%)
DT¼ 25Dt (%)
DT¼ 50Dt (%)
2 4 8 10 16
60.0 40.0 30.0 28.0 25.0
55.0 32.5 21.3 19.0 15.6
52.0 28.0 16.0 13.6 10.0
51.0 26.5 14.3 11.8 8.1
Notes: Np is the number of processors available for distribution; DT¼ Nfine Dt means that based on the two-level time-domain decomposition, each coarser time interval at the top level is divided into Nfine finer time intervals at the bottom level.
Table 5 Experimental efficiencies of the distributed ILPK algorithm to solve Example 2. Np
DT¼5Dt (%)
DT¼ 10Dt (%)
DT¼ 25Dt (%)
DT¼ 50Dt (%)
2 4 8 10 16
60.2 40.3 32.1 29.4 26.6
55.8 32.5 23.8 21.5 17.6
53.4 28.0 18.8 16.7 12.9
52.9 26.5 16.7 14.7 11.3
Fig. 3. Discrepancies of solutions of linear model obtained using the distributed ILPK algorithm with different parameters (for Example 1), notes: m is the number of parametric equations used in the numerical inverse Laplace transformation, Dt is the temporal mesh size.
Notes: Np is the number of processors available for distribution; DT¼ Nfine Dt means that based on the two-level time-domain decomposition, each coarser time interval at the top level is divided into Nfine finer time intervals at the bottom level.
LenovoTM PCs with PentiumTM 4 2.93 GHz processor, interconnected via Ethernet and CISCOTM 2950 switch. Experimental efficiency is evaluated using Eq. (15), with time costing of the distributed algorithm TC(parallel) divided by that of sequential algorithm TC(sequential). The sequential ILPK algorithm was run in the same computational environment with the same parameters. Experimental efficiency ¼
TCðparallelÞ TCðsequentialÞ
ð15Þ
Numerical data for the experimental efficiency is displayed in Table 5. Compare with the theoretical efficiency, the experimental efficiency is closer to but slightly greater than the theoretical value. This is due to the attribution of minor communication overheads of processors when solving the problem. 5.4. Parameter settings of the distributed ILPK The parameters involved in the distributed ILPK algorithm are the number of parametric equations m used in the numerical inverse Laplace transformation, the finer temporal mesh size Dt, and the coarser mesh size DT¼ Nfine Dt. It was reported that m should be set as an even number within 8, 10, y, 32 [14]. In this section, the accuracy of the algorithm due to these parameters is examined. To investigate the affects of these parameter settings on the accuracy of the algorithm for linear problems, Example 1 was solved by using the distributed ILPK algorithm with variations of m¼ 6,8,10,12,14,16,18, and Dt ¼0.1,0.01,0.001,0.0001. The total discrepancies of solutions obtained by the distributed algorithm and the analytical solutions were calculated and illustrated in Fig. 3. It is clear that, the discrepancy goes off smoothly when m Z10. In addition, Dt has little effect on discrepancies when it is no more than 0.01. To investigate the effect of these parameters of the distributed ILPK algorithm for nonlinear problems, the sum of the discrepancies between sequential ILPK for Example 2 are recorded in Fig. 4 with m¼6,8,10,12,14,16,18 and DT¼5Dt,10Dt,25Dt,50Dt, respectively. Results indicate that the value of m should be no less than 8. In addition, the larger DT leads to larger deviations to the results of the distributed ILPK algorithm.
Fig. 4. Discrepancies of solutions of nonlinear model obtained by distributed ILPK algorithm with different parameters (for Example 2), notes: m is the number of parametric equations used in the numerical inverse Laplace transformation, Dt is the finer temporal mesh size, and DT ¼ Nfine Dt is the coarser mesh size.
6. Conclusions In this paper, a time-domain decomposition methods are adapted for nonlinear pharmacokinetic models. The ILPK algorithm was implemented to solve pharmacokinetic models with an iterative inverse Laplace method in each time interval. A distributed ILPK algorithm based on the two-level time-domain decomposition is proposed to improve its efficiency. Solutions at top level on the coarser temporal mesh are obtained one by one, and then those at bottom level on the finer time intervals are calculated concurrently with initial solutions that have been obtained at the top level decomposition at suitable coarse time intervals. Experimental results show that the ILPK algorithm outperform the RK4 method for solving linear pharmacokinetic models and is able to obtain comparable solutions to analytical solutions. When solving nonlinear pharmacokinetic models, the discrepancies between solutions of the ILPK and the RK4 algorithms are very closer. It shows that the ILPK algorithm is a good candidate for both linear and nonlinear pharmacokinetic models. In the second part of this paper, a two-level decomposition of the time-domain to the distributed computing environments is presented. Experimental results show that the distributed algorithm is able to reduce computing time with regard to the use of more processors and larger time intervals DT¼Nfine Dt. Meanwhile, the coarser time interval DT affects the accuracy of the distributed ILPK algorithm. The authors recommend the relation between coarser
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
mesh and finer mesh as DT¼5Dt and 10Dt when the distributed ILPK algorithm produced solutions similar to the sequential ILPK algorithm. In conclusion, our experiments exploited accuracy–efficiency trade-off of the distributed algorithms. In addition, the proposed distributed algorithm is particularly suitable to problems that do not require detailed knowledge of each intermediate time when solving the pharmacokinetic models.
Conflict of interest statement No financial support from other organizations unless the cited in the Acknowledgments section, nor personal relationships, nor people, nor organisms biased the development of this research work.
227
[10] J.L. Fidalgo, W.K. Wong, Design issues for the Michaelis–Menten model, Journal of Theoretical Biology 215 (2002) 1–11. [11] J.C. Butcher, G. Wanner, Runge–Kutta methods: some historical notes, Applied Numerical Mathematics 22 (1-3) (1996) 113–151. [12] D.P. Gaver, Observing stochastic processes, and approximate transform inversion, Operations Research 14 (3) (1966) 444–459. [13] H. Stehfest, Numerical inversion of Laplace transforms, Communications of the ACM 13 (1970) 47–49. [14] B. Davies, B. Martin, Numerical inversion of the Laplace transform: a survey and comparison of methods, Journal of Computational Physics 33 (1) (1979) 1–32. [15] F. Yang, N. Li, A. Liu, Application of non-linear least square method and BP neural networks in the pharmacokinetics simulation of extravascular administration, Journal of Mathematical Medicine 20 (2) (2007) 200–202 in Chinese. [16] Y.F. Su, L.Y. Du, Plasma concentration of drugs obeying the Michaelis–Menten clearance kinetics in the case of extravascular administration, Herald of Medicine 25 (14) (2006) 357–359 in Chinese. [17] D. Barbolosi, A. Iliadis, Optimizing drug regimens in cancer chemotherapy: a simulation study using a PK-PD model, Computers in Biology and Medicine 31 (2001) 157–172.
Acknowledgments This research is supported by the Fundamental Research Funds for the Central Universities No. JUSRP10928 and innovation team project of Jiangnan University No. JNIRT0702. We would like to express our great appreciation to reviewers for their constructive suggestions and comments. References [1] K. Yamaoka, Y. Takakura, Analysis methods and recent advances in nonlinear pharmacokinetics from in vitro through in loci to in vivo, Drug Metabolism and Pharmacokinetics 19 (6) (2004) 397–406. [2] X.J. Huang, H.T. Xie, R.Y. Sun, Laplace transform and its applications in the pharmacokinetics, Chinese Journal of Clinical Pharmacology and Therapeutics 6 (1) (2001) 59–63 in Chinese. [3] J.H.E. Cartwright, O. Piro, The dynamics of Runge–Kutta methods, International Journal of Bifurcation and Chaos 2 (1992) 427–449. [4] E. Deeba, J.M.A. Yoon, Decomposition method for solving nonlinear systems of compartment models, Journal of Mathematical Analysis and Applications 266 (2002) 227–236. [5] J.E. Lagnese, G. Leugering, Time-domain decomposition of optimal control problems for the wave equation, Systems and Control Letters 48 (2003) 229–242. [6] P.C.A. De Haas, P.J. Zandbergen, The application of domain decomposition to time-domain computations of nonlinear water waves with a panel method, Journal of Computational Physics 129 (1996) 332–344. [7] C.H. Lai, A.K. Parrott, S.A. Rout, A distributed algorithm for European options with nonlinear volatility, Computers and Mathematics with Applications 49 (2005) 885–894. [8] G.J. Wang, Pharmacokinetics, Chemical Industry Press, 2005, pp. 88–120 in Chinese. [9] L. Rajaram, C.P. Tsokos, Statistical modeling of a pharmacokinetic system, Stochastic Analysis and Applications 24 (2006) 1061–1081.
Li LIU is an associate professor of School of IoT Engineering, Jiangnan University in China. She graduated with a B.Sc. degree in Computer Technology and its Applications and obtained a Ph.D. in the area of intelligent computing with applications in pharmacokinetics in 2008. She is an associate dean of Wuxi Fourth People’s Hospital, China. Her research interests include parallel computing and computational modeling for medical applications.
Choi-Hong LAI is Professor of Numerical Mathematics at the Department of Mathematical Sciences, University of Greenwich. He graduated with a first class B.Sc. degree in Mathematics and Engineering and obtained a Ph.D. in the area of numerical partial differential equations with applications in aerodynamics and parallel computing in 1985. He is the head of the research group Numerical and Applied Mathematics Unit within the School of Computing and Mathematical Sciences. His research interests include computational methods for medical applications, image processing and environmental applications.
Shao-Dan ZHOU is an associate chief pharmacist in the pharmaceutical preparation section of the Fourth People’s Hospital of Wuxi, Jiangsu, China, since 2008. She graduated with a bachelor’s degree in pharmacy from China Pharmaceutical University in 2000. Her present research is in the clinical trial.
Fen XIE is a pharmacist-in-charge in the pharmaceutical preparation section of the Fourth People’s Hospital of Wuxi, Jiangsu, China, since 2010. She graduated with a master’s degree in pharmacology from the Huaxi medical college of Sichuan University in 2006. Her present research is in the clinical pharmacy and bacterial resistance.
Rui LU is a postgraduate student of School of IoT Engineering, Jiangnan University, China. His current research interest is computational modeling for medical applications.
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Two-level time-domain decomposition based distributed method for numerical solutions of pharmacokinetic models Li Liu a,c,n, Choi-Hong Lai b, Shao-Dan Zhou c, Fen Xie c, Rui Lua a
School of IoT Engineering, Jiangnan University, Wuxi 214122, Jiangsu Province, China School of Computing and Mathematical Sciences, University of Greenwich, Old Royal Naval College, Greenwich, London SE10 9LS, UK c Suzhou University Affiliated Fourth People’s Hospital, Wuxi Fourth People’s Hospital, Wuxi 214062, Jiangsu Province, China b
a r t i c l e i n f o
abstract
Article history: Received 29 March 2010 Accepted 4 February 2011
In order to predict variations of drug concentration during a given period of time, numerical solutions of pharmacokinetic models need to be obtained efficiently. Analytical solutions of linear pharmacokinetic models are usually obtained using the Laplace transform and inverse Laplace tables. The derivations of solutions to complex nonlinear models are tedious, and such solution process may be difficult to implement as a robust software code. For nonlinear models, the fourth-order Runge–Kutta (RK4) is the most classical numerical method in obtaining approximate numerical solutions, which is impossible to be implemented in distributed computing environments without much modification. The reason is that numerical solutions obtained by using RK4 can only be computed in sequential time steps. In this paper, time-domain decomposition methods are adapted for nonlinear pharmacokinetic models. The numerical Inverse Laplace method for PharmacoKinetic models (ILPK) is implemented to solve pharmacokinetic models with iterative inverse Laplace transform in each time interval. The distributed ILPK algorithm, which is based on a two-level time-domain decomposition concept, is proposed to improve its efficiency. Solutions on the coarser temporal mesh at the top level are obtained one by one, and then those on the finer temporal mesh at the bottom level are calculated concurrently by using those initial solutions that have been obtained at the top level decomposition. Accuracy and efficiency of the proposed algorithm and its distributed equivalent are investigated by using several test models. Results indicate that the ILPK algorithm and its distributed equivalent are good candidates for both linear and nonlinear pharmacokinetic models. Crown Copyright & 2011 Published by Elsevier Ltd. All rights reserved.
Keywords: Pharmacokinetics Distributed computing Numerical inverse Laplace Runge-Kutta
1. Introduction Pharmacokinetic models are used to describe the absorption, distribution, metabolization and elimination (ADME) of drugs in vitro, in loci and in vivo [1]. In order to predict variations of drug concentrations during a given period of time, numerical solutions of the models need to be obtained efficiently. For linear models, the Laplace transform method supplemented with a lookup inverse Laplace table is often used to obtain the corresponding analytical solutions [2]. For nonlinear models, the fourth-order Runge–Kutta (RK4) [3] is a typical numerical method widely used to obtain approximate solutions, while the Adomian decomposition method has been reported to be used to obtain semianalytical solutions [4]. In practice, derivations of analytical and semi-analytical solutions are tedious and such solution process is
n Corresponding author at: School of IoT Engineering, Jiangnan University, Wuxi 214122, Jiangsu Province, China. E-mail address: [email protected] (L. Liu).
also difficult to be implemented as a robust software code. On the other hand, numerical solutions along the temporal axis calculated by using the RK4 method can only be obtained in sequential time steps. Inevitably, computing time of such temporal integration method becomes significant for nonlinear problems. The goal of this paper is to present an alternative approach to obtain numerical solutions of both linear and nonlinear pharmacokinetic models. At the same time, the efficiency of the method implemented using the distribution strategy is discussed and tested. Our approach is based on a time-domain decomposition method. There are a number of time-domain decomposition methods (TDDM). The key point of the time-domain decomposition is to split up the temporal axis into smaller intervals, which are used for the transmission of data across the intervals [5]. In this paper, a time-domain decomposition method is applied to solve a first-order nonlinear system. The ILPK algorithm is implemented to deal with the iterative inverse Laplace method in each time interval. One property that domain decomposition methods have is the suitability for parallel computing [6]. In our previous work, a two-level time-domain decomposition method
0010-4825/$ - see front matter Crown Copyright & 2011 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.compbiomed.2011.02.003
222
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
was applied to obtain numerical solutions of time-dependent nonlinear problems for European options [7]. In this paper, it is used to decouple the numerical calculation process of pharmacokinetic models. At the top level, the entire time-domain is divided into coarser time intervals. At the bottom level, each coarser time interval is further sub-divided into finer time intervals. Solutions on the coarser temporal mesh are obtained one by one in a sequential manner, and those on the finer temporal mesh are calculated concurrently with initial solutions that have been obtained at the top level decomposition. Further distributed properties are also induced into the algorithm at the coarser temporal computation by using Laplace transformation. In the numerical experiments, the sequential and the distributed methods are implemented and tested for several models chosen from pharmacokinetics and both the accuracy and the efficiency of the proposed methods are investigated. Parameter settings of these methods are also discussed.
problem of solving linear or nonlinear model is converted into obtaining solutions of Eq. (3). In the next section, numerical methods are given in details.
3. Numerical methods to solve pharmacokinetic models 3.1. The Runge–Kutta methods The Runge–Kutta methods represent a class of numerical methods in obtaining approximate numerical solutions of differential equations related to initial value problems [11]. The methods rely on the Taylor series expansion but do not require the computation of high order derivatives of the right hand side function of the differential equations. The fourth-order Runge– Kutta (RK4) is a typical numerical method widely used to obtain approximate solutions, which requires four evaluations of updates per time step. Its error per time step Dt is of the order of Dt5, while the total accumulated error is of the order of Dt4.
2. The pharmacokinetic models Compartment model [8] is a classical type of pharmacokinetic models and views the body as a system composed of many compartments taking into variations of drug concentration. The distribution of a drug is described usually as transportation of drug concentration between compartments and is represented mathematically by a set of differential equations. For example, an n-compartment model in Ref. [9] is defined as follows: x_ i ðtÞ ¼
n X
kij xi ðtÞ þ
j¼1 iaj
n X
kji xj ðtÞ,
i,j ¼ 1,2,3,. . .,n
ð1Þ
j¼1 i aj
with initial conditions xi(0)Z0. Here xi(t) represents the drug concentration in the ith compartment at time t, kij is the proportional constant associated with the rate of absorption or elimination from the ith compartment to the jth compartment. In addition, there are many drugs which are eliminated by enzymes and transporters expressed mainly in the intestine, liver and kidney, and are bound to proteins in the blood and tissues. Since the numbers of enzymes, transporters and binding proteins are limited, the disposition of drugs is capacity-limited. In this context, the Michaelis–Menten model is often used to describe the rapid equilibrium kinetics of enzyme systems [10]. It is also embedded into compartment models to represent the rate of change from one compartment to another. Here is an example of one compartment model with Michaelis–Menten elimination. The drug is given by intravenous injection. Let Vm be the maximum elimination rate, and Km be the Michaelis–Menten constant. The coefficient Vm/(Km + x(t)) represents the capacity-limited clearance. The variation of drug concentration described by Michaelis–Menten equation is defined as below _ ¼ xðtÞ
Vm xðtÞ Km þ xðtÞ
ð2Þ
A generalized mathematical description of pharmacokinetic models capable of describing both linear and nonlinear compartment models can be written as X_ ¼ AX þHðXÞ þ FðtÞ, t A ½0,T
ð3Þ
where X is the concentration matrix of compartments and X_ ¼ dX=dt. A multi-compartment model composes three parts: the linear part AX represents first-order kinetic distribution between compartments where A is a coefficient matrix, the nonlinear part H(X) representing nonlinear absorption or elimination and the external input function F(t) which is used to express the administration protocol of drugs. As a result, the
3.2. Numerical inverse Laplace and time-domain decomposition The Laplace transform of the time dependent function x(t) is an integral transform defined by the integral Z 1 elt xðtÞdt ¼ UðlÞ ð4Þ LðxðtÞÞ 0
It is particularly useful in solving linear ordinary differential equations where an ordinary differential equation with temporal derivatives is transformed into a form that does not involve time variable. In this paper the authors examine an application of Laplace transformation to nonlinear differential equations. There are many different methods for numerically inverting the Laplace transform. Among them, the Graver–Stehfest inverse Laplace method which was proposed by Graver [12] and derived by Stehfest [13] gave good accuracy on a fairly wide range of problems [14]. It should be noted that the Laplace integral is known to exist on the condition that xðtÞ be piecewise continuous on every finite interval in t Z 0 and satisfy xðtÞ r Meat for certain constants M and a. It is clear that drug concentrations are under this condition. Hence, taking Laplace transform of Eq. (3) leads to ð5Þ UðlÞ ¼ ½lIA1 ½LðHðXÞÞ þLðFðtÞÞ þ Xðt0 Þ subject to the constraint lIA a0. Here l A flj ,j ¼ 1,2,. . .,mg is a finite set of transformation parameters defined as below
lj ¼ j
ln 2 , T
j ¼ 1,2,. . .,m
ð6Þ
In order to retrieve X, the numerical inverse Laplace transform due to Graver–Stehfest method as listed in Eq. (7) may be used X
m ln 2 X U w UUðlj Þ T j¼1 j
ð7Þ
where the weighting factor wj is given by 2 3 minðj,m=2Þ X im=2 ð2iÞ! j þ ðm=2Þ 4 5 U wj ¼ ð1Þ ðm=2iÞ!i!ði1Þ!ðjiÞ!ð2ijÞ! i ¼ ð1 þ jÞ=2 In the Graver–Stehfest method, m is the only parameter to be set. As defined in Eq. (7), the approximation of Laplace transformed X is obtained by integrating m different parametric functions of Uðlj Þ. The higher the value of m, the more computing time is required to calculate X but not necessarily more accurate. Note that it is possible to compute these parametric functions simultaneously using a distributed computing environment.
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
The reason is that these parametric functions are independent from each other. To solve the nonlinear part of the model, time-domain decomposition is applied firstly. The key point of the time-domain decomposition is to split up the temporal axis ½0,T into smaller intervals Ip ¼ ½tp ,tp þ 1 with t0 ¼ 0, tN ¼ T, tp otp þ 1 , and p ¼ 0,1,. . .,N1, in order to transmit data across the boundary of Ip . As such the numerical inverse Laplace method has to be done sequentially in each time interval. As in Eq. (8), HðXÞ is computed using the linearization X~ , which is updated within an iterative process until convergence is achieved for HðXÞ within that time interval. Each step of the iterative update process involves a numerical solution to the equation X_ ðtÞ ¼ AX þ HðX~ Þ þ FðtÞ,
t A ½tp ,tp þ 1
ð8Þ
Let Xðtp Þ and Xðtp þ 1 Þ be the numerical solutions of Eq. (3) at tp and tp þ 1 . A Laplace transform is applied to Eq. (8), being defined on t A ½tp ,tp þ 1 , in its differential form which leads to UðX~ , lÞ ¼ ½lIA1 ½HðX~ Þ=l þ LðFðtÞÞ þXðtp Þ. In this paper, the iterative inverse Laplace to obtain numerical solution Xðtp þ 1 Þ of nonlinear pharmacokinetic models, using Xðtp Þ as the initial approximation to X~ is implemented. The approximation X~ of Xðtp þ 1 Þ is updated with X ðkÞ ðtp þ 1 Þ at the kth step iteration until the convergence criteria :X ðkÞ ðtp þ 1 ÞX ðk1Þ ðtp þ 1 Þ: o e is achieved. Here is the main function the ILPK algorithm used to compute the solution X. FUNCTION X(tp + 1)¼ILPK(X(tp)) Initial approximation: X(0)(tp + 1) ¼X(tp); k¼0; // Iterative Inverse Laplace in Ip ¼[tp,tp + 1] ITERATE X~ ¼ X ðkÞ ðtp þ 1 Þ; Compute simultaneous for j¼ 1, 2, y, m: UðX~ , lj Þ ¼ ½lj IA1 ½HðX~ Þ=lj þ LðFðtÞÞ þ Xðtp Þ, where lj ¼ j tp þln1 2tp ; Compute X ðk þ 1Þ ðtp þ 1 Þ
ln 2 tp þ 1 tp
" j þ ðm=2Þ
where wj ¼ ð1Þ
U
m P
wj UðX~ , lj Þ,
j¼1
minðj,m=2Þ P i ¼ ð1 þ jÞ=2
im=2 ð2iÞ! ðm=2iÞ!i!ði1Þ!ðjiÞ!ð2ijÞ!
# ;
k ¼k+1; UNTIL :X ðkÞ ðtp þ 1 ÞX ðk1Þ ðtp þ 1 Þ:o e; RETURN X(k)(tp + 1); END-FUNCTION
If H(X) ¼0, Eq. (3) represents a linear model, and its solution at time t ¼tp + 1, p ¼0, 1, ..., N 1 can be obtained by X(tp + 1)¼ ILPK(X(t0)) simultaneously when distributed to N processors. However, for the nonlinear model, solutions can only be obtained sequentially since computing X(tp + 1) is very depended on its initial solution X(tp). Section 4 proposes a method of improving the efficiency of the ILPK algorithm under distributed environments—two-level time-domain decomposition.
4. Distributed algorithms based on two-level time-domain decomposition In our previous work, a two-level time-domain decomposition method was used to decouple time-dependent nonlinear problems for European options. The main idea of the two-level timedomain decomposition is to divide the entire time-domain [0,T] into smaller time intervals. At the bottom level, [0,T] is equally divided into finer time intervals Ibp ¼[tp,tp + 1], with, tp + 1 ¼tp + Dt, p ¼0, 1, ..., N 1, t0 ¼0, tN ¼T. At the top level, Nfine pieces of finer
223
Fig. 1. Illustration of the process of the distributed ILPK based on the two-level time-domain decomposition with DT¼ 5Dt
time intervals are integrated into a coarse time interval, Itq ¼ [Tq,Tq + 1], with Tq + 1 ¼Tq +Nfine Dt, q¼0, 1, ..., M 1, M ¼ N=Nfine and T0 ¼t0. Based on the two-level time-domain decomposition, X(Tq), q¼0, 1, ..., M 1 on the coarser temporal mesh are computed one by one in a sequential manner, and are being used as initial conditions of solutions on the finer time intervals. So remaining solutions on the finer intervals X(Tq + ti),ti ¼i Dt, and those X(Tk + ti), k aq, i¼ 1, ..., Nfine 1 can be obtained concurrently. In short, sequential X(tp) is mapped into two-level time domain X(Tq + ti), q ¼ p=Nfine , i ¼ kmodNfine . Given DT¼5Dt, the process of the distributed ILPK based on the two-level timedomain decomposition is illustrated in Fig. 1. For example, after X(t5), X(t10), X(t15) have been calculated sequentially, X(t1) is obtained with initial condition X(t0) followed by X(t2), X(t3), X(t4) and parallelized with X(t6), X(t11)y The pseudo-code of the distributed ILPK algorithm is described as below. FUNCTION Distributed ILPK IF H(X)¼0, // for linear models, do one-level time-domain decomposition THEN Distribute p ¼1, y, N tp to processors, X(tp) ¼ILPK(X(t0)); END-Distribute and receive X(tp). ELSE // for nonlinear models, do two-level time-domain decomposition T0 ¼t0; DT¼Nfine Dt; M ¼ N=Nfine ; FOR q¼0 To M 1, do coarser temporal mesh on tA[Tq,Tq + 1] sequentially Tq + 1 ¼Tq + DT; // top level time interval X(Tq + 1) ¼ILPK(X(Tq)); END-FOR Distribute q¼0 To M 1,X(Tq) to processors, For i¼1 To Nfine 1, do finer temporal mesh on tA[Tq + ti 1,Tq + ti] ti ¼i Dt; // bottom level time interval X(Tq + ti)¼ ILPK(X(Tq + ti 1)); END-FOR END-Distribute and receive X(Tq + t1), y, X(Tq + tNfine 1 ). END-IF END-FUNCTION
From the pseudo-code above we can see that if we replace the ILPK algorithm with the RK4, the computing of solutions can also be distributed based on the two-level time-domain decomposition method. In the next section, numerical experiments are provided to examine the accuracy and the efficiency of the distributed ILPK and RK4 algorithms.
5. Experiments and discussion In our experiments, the ILPK and the RK4 methods are programmed with MatlabTM software. In Section 5.1, linear,
224
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
nonlinear pharmacokinetic models, as well as an integrated pharmacokinetic–pharmacodynamic model are solved by using the ILPK algorithm. In Section 5.2, the accuracy of the distributed algorithms with variation of two-level time-domain distributed strategy DT ¼Nfine Dt is investigated. Theoretical and experimental efficiencies are calculated and observed with Matlab distributed computing platform in Section 5.3. Finally, parameter setting of the distributed algorithms is discussed in Section 5.4.
Solutions X(ILPK) and X(RK4) in tA[0,0.05] and their discrepancies are shown in Table 2. The discrepancies :XðILPKÞXðRK4Þ: are less than 4E-6.
5.1. Accuracy analysis of the ILPK algorithm
A more complicated model was a pharmacokinetic– pharmacodynamic (PK–PD) model, which was built by Barbolosi [17] to predict cancer chemotherapy. Here, we extend the application of the ILPK method to solve this model. c1(t) is the concentration of administered anticancer drug in the plasma. The mass balance for the two-compartment model leads to 8 1 ðtÞ < dcdt ¼ ðkE1 þ k12 ÞUc1 ðtÞ þ uðtÞ V1 ð11Þ : dc2 ðtÞ ¼ k12 U V1 Uc1 ðtÞkE2 Uc2 ðtÞ V2 dt
Example 1. A linear model for extra vascular administration A linear pharmacokinetic model for oral administration proposed in literature [15] is rewritten as Eq. (3) and displayed in Eq. (9). Here, H(X)¼0, F(t)¼0. Given initial dose, drug went into " # x1 compartments with first-order absorption process. X ¼ x2 represented the concentrations of the two compartments. The rate constants of the differential equation were k12 ¼2.008 and k21 ¼0.1855. The initial drug concentrations in the two compartments at t ¼0 were x1(0) ¼103.8613 and x2(0) ¼0. " # k12 0 _ X¼ X ð9Þ k21 k12 The classical method to obtain analytical solutions x1, x2 of the linear compartment model is the Laplace transform supplemented with a look-up inverse Laplace table. Here, analytical solutions of x1, x2 were deduced. The parameter of both the ILPK and RK4 methods were set as m¼ 16, and time step Dt ¼0.001. Approximate solutions obtained by the ILPK and RK4 method were denoted as x1(ILPK), x2(ILPK), x1(RK4) and x2(RK4). Discrepancies between analytical solutions and those obtained by the ILPK and the RK4 method :x1 x1 ðRK4Þ:, :x2 x2 ðRK4Þ:, :x1 x1 ðILPKÞ: and :x2 x2 ðILPKÞ: are recorded in Table 1. It is clear that the approximate solutions obtained by the ILPK algorithm are closer to analytical solutions than the RK4 method. Example 2. Nonlinear elimination
pharmacokinetic
model
for
drug
A model obeying the first order and Michaelis–Menten elimination kinetics in the case of extra vascular administration was displayed in Eq. (10), which was solved by the RK4 method [16]. Here, the model with constants ka ¼0.5 h 1, Vd/F¼40 L, 1 D ¼150 mg, tlag ¼0, Vm ¼ 0:8 mg L1 h , Km ¼10 mg L 1 is solved by using the ILPK method. In the numerical experiments the parameters of the ILPK algorithm were set as m¼10, Dt ¼0.0001.
dX Vm X FD ka ðttlag Þ ¼ kX þ ka e dt Km þX Vd
ð10Þ
Example 3. Pharmacokinetic–pharmacodynamic model of cancer chemotherapy
where 8 > < 0:205 uðtÞ ¼ 0:1 > : 0
for
0 rt o 1
for
1r t o22
for
22 rt
is the drug delivery by intravenous infusion. Assuming that the growth of tumor behaved according to a Gompertz-type growth equation, n(t) denotes the number of tumor cells at time t, and H(.) is a Heaviside step function dnðtÞ ¼ lUnðtÞUln y=nðtÞ kU½c2 ðtÞCMIN UnðtÞUHðc2 ðtÞCMIN Þ dt
ð12Þ
Table 2 Solutions of Example 2 obtained using the ILPK and RK4 methods. t(h)
X(ILPK)
X(RK4)
:XðILPKÞXðRK4Þ:
0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050
0.0093613 0.0186954 0.0280025 0.0372826 0.0465358 0.0557622 0.0649618 0.0741347 0.0832810 0.0924008
0.0093610 0.0186960 0.0280030 0.0372830 0.0465370 0.0557630 0.0649630 0.0741360 0.0832830 0.0924030
3.53E-7 7.05E-7 1.06E-6 1.41E-6 1.75E-6 2.10E-6 2.45E-6 2.79E-6 3.14E-6 3.48E-6
Notes: X(ILPK), X(RK4) were approximate solutions obtained using the ILPK and RK4 methods, respectively. :XðILPKÞXðRK4Þ: were discrepancies between these solutions.
Table 1 Discrepancies between analytical solutions and those obtained by the ILPK and RK4 methods. t(h)
x1
x2
:x1 x1 ðRK4Þ:
:x2 x2 ðRK4Þ:
:x1 x1 ðILPKÞ:
:x2 x2 ðILPKÞ:
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
84.7960 69.3696 56.7496 46.4255 37.9796 31.0702 25.4178 20.7937 17.0108 13.9161
18.8818 33.8143 45.6926 55.0786 62.4318 68.1279 72.4744 75.7224 78.0775 79.7078
0.1704 0.1394 0.1141 0.0933 0.0763 0.0625 0.0511 0.0418 0.0342 0.0280
0.1670 0.1332 0.1056 0.0831 0.0648 0.0498 0.0376 0.0278 0.0197 0.0132
3E-06 2E-06 1E-05 4E-06 1.1E-05 1.9E-05 2E-06 4.2E-05 8.8E-05 1.7E-04
0 1E-06 1E-06 4E-06 9E-06 1.2E-05 0 3.8E-05 1.1E-04 1.8E-04
Notes: x1, x2 were analytical solutions of Example 1. :x1 x1 ðRK4Þ:, :x2 x2 ðRK4Þ:, :x1 x1 ðILPKÞ:, :x2 x2 ðILPKÞ: were discrepancies between analytical solutions and those obtained by the ILPK and RK4 methods.
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
225
In addition, dynamical toxicities induced by drug exposures is also taken into account in a cumulative process expressing myelosuppression in the bone marrow site, which is formulated as dwðtÞ ¼ rc vUwðtÞmUwðtÞUc1 ðttÞ dt
ð13Þ
In the numerical experiments the constants of the model were V1 ¼ 20 L, kE1 ¼ 0.48 day 1, k12 ¼0.02 day 1, V2 ¼3 L, kE2 ¼0.2 day 1, l ¼3 10 3 day 1, y ¼1012, k ¼ 10 g1 L day1 , CMIN ¼1 mg mL 1, rc ¼0.8 10 3 L 1 day 1, v¼0.1 day 1, m ¼50 g 1 L day 1, t ¼10 days. The initial values of the model were c1(0)¼0 g L 1, c2(0)¼ 0 g L 1, n(0)¼n0, w(0)¼w0, no ¼5 1010 cells, wo ¼rc/v, tA [0,45] days. To present the variations of drug concentrations, the number of tumor cells, and the number of white blood cells during 45 days, the ILPK and the RK4 method were used. Since there was time delay for toxicities, Eq. (13) was solved after Eqs. (11) and (12) had been done. The parameters of the ILPK algorithm were set as m¼10, and Dt ¼0.005. Discrepancies between approximate solutions obtained by the ILPK and RK4 methods are recorded in Table 3. It indicates that solutions obtained by the ILPK algorithm are very closer to those of the RK4 method.
Fig. 2. Comparison of discrepancies between solutions of Example 3 obtained using the distributed ILPK, RK4 algorithms and using their sequential ones.
those of the distributed ILPK algorithm. It also indicates that the accuracy of the distributed ILPK algorithm is little affected when DT¼5Dt and 10Dt. 5.3. Efficiency analysis of the distributed algorithms
5.2. Accuracy analysis of the distributed algorithms It is clear that the time interval Dt¼tp + 1 tp is one main contributing factor to the accuracy of solutions obtained by both the sequential ILPK and the RK4 methods, since solution X(tp + 1) is obtained by using X(tp) as the initial approximation. Compared to its sequential equivalent, the distributed ILPK and RK4 methods calculate solutions in two phases. Firstly, solutions at the bottom level (finer intervals) are dependent on that at the top level (coarse intervals) solutions. It can be imagined that with a larger coarse interval DT¼Nfine Dt, discrepancies between solutions of the distributed algorithm and its sequential equivalent become greater. In this paper we investigate how far the accuracy of the distributed algorithms is being affected by the distributed strategy represented as coarse interval DT¼ Nfine Dt. The distributed ILPK and the distributed RK4 algorithms were run to solve Example 3 with DT¼ 5Dt,10Dt,25Dt,50Dt, respectively. The parameters of the algorithms used in the numerical experiments were m¼10 and Dt¼0.01. Accumulated discrepancies for 45 days between solutions obtained by using the distributed and its sequential algorithm are illustrated in Fig. 2. When DT increases from 5Dt to 50Dt, discrepancies between solutions obtained by the distributed algorithms and their sequential ones goes up quickly. In addition, the discrepancies between solutions of the distributed RK4 and its sequential version increased a little more quickly than Table 3 Solutions of Example 3 obtained by using the ILPK and RK4 methods. t(day)
c1(ILPK)
c2(ILPK)
n(ILPK)
w(ILPK)
d(c1)
d(c2)
d(n)
d(w)
5 10 15 20 25 30 35 40 45
0.00970 0.00998 0.01000 0.01000 0.00218 0.00018 0.00001 0.00000 0.00000
0.00361 0.00550 0.00624 0.00651 0.00501 0.00212 0.00080 0.00030 0.00011
9.52 7.92 6.22 4.80 3.75 3.38 3.37 3.43 3.48
0.08000 0.08000 0.01945 0.01370 0.01336 0.01333 0.02094 0.04034 0.05545
7E-05 2E-06 4E-06 3E-06 1E-04 1E-05 1E-05 0 0
7E-05 4E-05 1E-05 3E-06 8E-05 4E-05 2E-05 0 0
0.02 0.03 0.03 0.03 0.03 0.04 0.04 0.05 0.04
2E-06 1E-06 3E-04 2E-05 1E-05 1E-06 4E-04 4E-04 3E-04
Notes: c1(ILPK), c2(ILPK), n(ILPK), w(ILPK) were approximate solutions of Example 3 obtained by the ILPK method. dðc1 Þ ¼ :c1 ðILPKÞc1 ðRK4Þ:, dðc2 Þ ¼ :c2 ðILPKÞ c2 ðRK4Þ:, dðnÞ ¼ :nðILPKÞnðRK4Þ:, dðwÞ ¼ :wðILPKÞwðRK4Þ: were discrepancies between approximate solutions of Example 3 obtained using the ILPK and RK4 methods.
The execution time of a distributed algorithm composes of computing time and communication overheads. The latter is governed by the computing network environments. When problems to be handled are more complex and can be parallelized without much communication between different problems, the communication overheads can be neglected. Defining one work unit as the computational work required for solving the pharmacokinetic model X(tp + 1)¼ ILPK(X(tp)) in a given time interval [tp,tp + 1]. The total parallel work unit of the distributed ILPK algorithm is composed of the total work unit at the top level and the distributed work unit at the bottom level executed in each processor. In this paper, the theoretical efficiency of the distributed algorithm is defined as Eq. (14), with the number of parallel work unit NWU(parallel) divided by the number of sequential work unit NWU(sequential) it required. The smaller values of efficiency show that the algorithm is more efficient. It gives us a measure to evaluate how much time saving is achieved when using the distributed strategy without the effect of communication overheads. Theoretical efficiency ¼
NWUðparallelÞ NWUðsequentialÞ
ð14Þ
For linear models, X(tp + 1)¼ILPK(X(tp)), p¼0, 1, ..., N 1, the efficiency can be calculated simultaneously in processors. Given the distributed algorithm is executed in Np processors, its theoretical efficiency to solve linear model is 1/Np. When solving nonlinear problems with two-level time-domain decomposition DT¼Nfine Dt, the NWU(parallel) is composed of n=Nfine work units of initial conditions and those on finer time intervals in each processor n (Nfine 1)/(NfineNp). So, the theoretical efficiency of the distributed ILPK algorithm for solving nonlinear problem is efficiency¼1/Nfine + 1/Np 1/(NfineNp), which is inversely proportional to the number of processors Np and the coefficient of coarser time interval Nfine. Table 4 shows the theoretical efficiency of the distributed ILPK calculated with a range of coarse time intervals DT¼5Dt, 10Dt,25Dt,50Dt, and with different number of processors Np ¼2,4,8,10,16 for the problem described in Example 2. We can see that the distributed algorithm saves much more computing time with regard to more processors and bigger time intervals. In this paper, experimental efficiency with various Nfine, Np was also observed. The distributed environment was set up with MatlabTM Distributed Computing Engine (MDCE). The performance analysis was obtained on a PC cluster consisting of 16
226
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
Table 4 Theoretical efficiencies of the distributed ILPK algorithm to solve Example 2. Np
DT¼5Dt (%)
DT¼ 10Dt (%)
DT¼ 25Dt (%)
DT¼ 50Dt (%)
2 4 8 10 16
60.0 40.0 30.0 28.0 25.0
55.0 32.5 21.3 19.0 15.6
52.0 28.0 16.0 13.6 10.0
51.0 26.5 14.3 11.8 8.1
Notes: Np is the number of processors available for distribution; DT¼ Nfine Dt means that based on the two-level time-domain decomposition, each coarser time interval at the top level is divided into Nfine finer time intervals at the bottom level.
Table 5 Experimental efficiencies of the distributed ILPK algorithm to solve Example 2. Np
DT¼5Dt (%)
DT¼ 10Dt (%)
DT¼ 25Dt (%)
DT¼ 50Dt (%)
2 4 8 10 16
60.2 40.3 32.1 29.4 26.6
55.8 32.5 23.8 21.5 17.6
53.4 28.0 18.8 16.7 12.9
52.9 26.5 16.7 14.7 11.3
Fig. 3. Discrepancies of solutions of linear model obtained using the distributed ILPK algorithm with different parameters (for Example 1), notes: m is the number of parametric equations used in the numerical inverse Laplace transformation, Dt is the temporal mesh size.
Notes: Np is the number of processors available for distribution; DT¼ Nfine Dt means that based on the two-level time-domain decomposition, each coarser time interval at the top level is divided into Nfine finer time intervals at the bottom level.
LenovoTM PCs with PentiumTM 4 2.93 GHz processor, interconnected via Ethernet and CISCOTM 2950 switch. Experimental efficiency is evaluated using Eq. (15), with time costing of the distributed algorithm TC(parallel) divided by that of sequential algorithm TC(sequential). The sequential ILPK algorithm was run in the same computational environment with the same parameters. Experimental efficiency ¼
TCðparallelÞ TCðsequentialÞ
ð15Þ
Numerical data for the experimental efficiency is displayed in Table 5. Compare with the theoretical efficiency, the experimental efficiency is closer to but slightly greater than the theoretical value. This is due to the attribution of minor communication overheads of processors when solving the problem. 5.4. Parameter settings of the distributed ILPK The parameters involved in the distributed ILPK algorithm are the number of parametric equations m used in the numerical inverse Laplace transformation, the finer temporal mesh size Dt, and the coarser mesh size DT¼ Nfine Dt. It was reported that m should be set as an even number within 8, 10, y, 32 [14]. In this section, the accuracy of the algorithm due to these parameters is examined. To investigate the affects of these parameter settings on the accuracy of the algorithm for linear problems, Example 1 was solved by using the distributed ILPK algorithm with variations of m¼ 6,8,10,12,14,16,18, and Dt ¼0.1,0.01,0.001,0.0001. The total discrepancies of solutions obtained by the distributed algorithm and the analytical solutions were calculated and illustrated in Fig. 3. It is clear that, the discrepancy goes off smoothly when m Z10. In addition, Dt has little effect on discrepancies when it is no more than 0.01. To investigate the effect of these parameters of the distributed ILPK algorithm for nonlinear problems, the sum of the discrepancies between sequential ILPK for Example 2 are recorded in Fig. 4 with m¼6,8,10,12,14,16,18 and DT¼5Dt,10Dt,25Dt,50Dt, respectively. Results indicate that the value of m should be no less than 8. In addition, the larger DT leads to larger deviations to the results of the distributed ILPK algorithm.
Fig. 4. Discrepancies of solutions of nonlinear model obtained by distributed ILPK algorithm with different parameters (for Example 2), notes: m is the number of parametric equations used in the numerical inverse Laplace transformation, Dt is the finer temporal mesh size, and DT ¼ Nfine Dt is the coarser mesh size.
6. Conclusions In this paper, a time-domain decomposition methods are adapted for nonlinear pharmacokinetic models. The ILPK algorithm was implemented to solve pharmacokinetic models with an iterative inverse Laplace method in each time interval. A distributed ILPK algorithm based on the two-level time-domain decomposition is proposed to improve its efficiency. Solutions at top level on the coarser temporal mesh are obtained one by one, and then those at bottom level on the finer time intervals are calculated concurrently with initial solutions that have been obtained at the top level decomposition at suitable coarse time intervals. Experimental results show that the ILPK algorithm outperform the RK4 method for solving linear pharmacokinetic models and is able to obtain comparable solutions to analytical solutions. When solving nonlinear pharmacokinetic models, the discrepancies between solutions of the ILPK and the RK4 algorithms are very closer. It shows that the ILPK algorithm is a good candidate for both linear and nonlinear pharmacokinetic models. In the second part of this paper, a two-level decomposition of the time-domain to the distributed computing environments is presented. Experimental results show that the distributed algorithm is able to reduce computing time with regard to the use of more processors and larger time intervals DT¼Nfine Dt. Meanwhile, the coarser time interval DT affects the accuracy of the distributed ILPK algorithm. The authors recommend the relation between coarser
L. Liu et al. / Computers in Biology and Medicine 41 (2011) 221–227
mesh and finer mesh as DT¼5Dt and 10Dt when the distributed ILPK algorithm produced solutions similar to the sequential ILPK algorithm. In conclusion, our experiments exploited accuracy–efficiency trade-off of the distributed algorithms. In addition, the proposed distributed algorithm is particularly suitable to problems that do not require detailed knowledge of each intermediate time when solving the pharmacokinetic models.
Conflict of interest statement No financial support from other organizations unless the cited in the Acknowledgments section, nor personal relationships, nor people, nor organisms biased the development of this research work.
227
[10] J.L. Fidalgo, W.K. Wong, Design issues for the Michaelis–Menten model, Journal of Theoretical Biology 215 (2002) 1–11. [11] J.C. Butcher, G. Wanner, Runge–Kutta methods: some historical notes, Applied Numerical Mathematics 22 (1-3) (1996) 113–151. [12] D.P. Gaver, Observing stochastic processes, and approximate transform inversion, Operations Research 14 (3) (1966) 444–459. [13] H. Stehfest, Numerical inversion of Laplace transforms, Communications of the ACM 13 (1970) 47–49. [14] B. Davies, B. Martin, Numerical inversion of the Laplace transform: a survey and comparison of methods, Journal of Computational Physics 33 (1) (1979) 1–32. [15] F. Yang, N. Li, A. Liu, Application of non-linear least square method and BP neural networks in the pharmacokinetics simulation of extravascular administration, Journal of Mathematical Medicine 20 (2) (2007) 200–202 in Chinese. [16] Y.F. Su, L.Y. Du, Plasma concentration of drugs obeying the Michaelis–Menten clearance kinetics in the case of extravascular administration, Herald of Medicine 25 (14) (2006) 357–359 in Chinese. [17] D. Barbolosi, A. Iliadis, Optimizing drug regimens in cancer chemotherapy: a simulation study using a PK-PD model, Computers in Biology and Medicine 31 (2001) 157–172.
Acknowledgments This research is supported by the Fundamental Research Funds for the Central Universities No. JUSRP10928 and innovation team project of Jiangnan University No. JNIRT0702. We would like to express our great appreciation to reviewers for their constructive suggestions and comments. References [1] K. Yamaoka, Y. Takakura, Analysis methods and recent advances in nonlinear pharmacokinetics from in vitro through in loci to in vivo, Drug Metabolism and Pharmacokinetics 19 (6) (2004) 397–406. [2] X.J. Huang, H.T. Xie, R.Y. Sun, Laplace transform and its applications in the pharmacokinetics, Chinese Journal of Clinical Pharmacology and Therapeutics 6 (1) (2001) 59–63 in Chinese. [3] J.H.E. Cartwright, O. Piro, The dynamics of Runge–Kutta methods, International Journal of Bifurcation and Chaos 2 (1992) 427–449. [4] E. Deeba, J.M.A. Yoon, Decomposition method for solving nonlinear systems of compartment models, Journal of Mathematical Analysis and Applications 266 (2002) 227–236. [5] J.E. Lagnese, G. Leugering, Time-domain decomposition of optimal control problems for the wave equation, Systems and Control Letters 48 (2003) 229–242. [6] P.C.A. De Haas, P.J. Zandbergen, The application of domain decomposition to time-domain computations of nonlinear water waves with a panel method, Journal of Computational Physics 129 (1996) 332–344. [7] C.H. Lai, A.K. Parrott, S.A. Rout, A distributed algorithm for European options with nonlinear volatility, Computers and Mathematics with Applications 49 (2005) 885–894. [8] G.J. Wang, Pharmacokinetics, Chemical Industry Press, 2005, pp. 88–120 in Chinese. [9] L. Rajaram, C.P. Tsokos, Statistical modeling of a pharmacokinetic system, Stochastic Analysis and Applications 24 (2006) 1061–1081.
Li LIU is an associate professor of School of IoT Engineering, Jiangnan University in China. She graduated with a B.Sc. degree in Computer Technology and its Applications and obtained a Ph.D. in the area of intelligent computing with applications in pharmacokinetics in 2008. She is an associate dean of Wuxi Fourth People’s Hospital, China. Her research interests include parallel computing and computational modeling for medical applications.
Choi-Hong LAI is Professor of Numerical Mathematics at the Department of Mathematical Sciences, University of Greenwich. He graduated with a first class B.Sc. degree in Mathematics and Engineering and obtained a Ph.D. in the area of numerical partial differential equations with applications in aerodynamics and parallel computing in 1985. He is the head of the research group Numerical and Applied Mathematics Unit within the School of Computing and Mathematical Sciences. His research interests include computational methods for medical applications, image processing and environmental applications.
Shao-Dan ZHOU is an associate chief pharmacist in the pharmaceutical preparation section of the Fourth People’s Hospital of Wuxi, Jiangsu, China, since 2008. She graduated with a bachelor’s degree in pharmacy from China Pharmaceutical University in 2000. Her present research is in the clinical trial.
Fen XIE is a pharmacist-in-charge in the pharmaceutical preparation section of the Fourth People’s Hospital of Wuxi, Jiangsu, China, since 2010. She graduated with a master’s degree in pharmacology from the Huaxi medical college of Sichuan University in 2006. Her present research is in the clinical pharmacy and bacterial resistance.
Rui LU is a postgraduate student of School of IoT Engineering, Jiangnan University, China. His current research interest is computational modeling for medical applications.