Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
Contents lists available at ScienceDirect
Journal of Pharmacological and Toxicological Methods journal homepage: www.elsevier.com/locate/jpharmtox
Original article
ADVAN-style analytical solutions for common pharmacokinetic models Ahmad Y. Abuhelwa, David J.R. Foster, Richard N. Upton ⁎ Australian Centre for Pharmacometrics and Sansom Institute, School of Pharmacy and Medical Sciences, University of South Australia, SA 5000, Australia
a r t i c l e
i n f o
Article history: Received 30 October 2014 Received in revised form 25 March 2015 Accepted 27 March 2015 Available online 2 April 2015 Keywords: ADVAN Analytical solutions Intravenous Methods Modelling NONMEM Oral Pharmacokinetics R-language Simulations
a b s t r a c t Introduction: The analytical solutions to compartmental pharmacokinetic models are well known, but have not been presented in a form that easily allows for complex dosing regimen and changes in covariate/parameter values that may occur at discrete times within and/or between dosing intervals. Methods: Laplace transforms were used to derive ADVAN-style analytical solutions for 1, 2, and 3 compartment pharmacokinetic linear models of intravenous and first-order absorption drug administration. The equations calculate the change in drug amounts in each compartment of the model over a time interval (t; t = t2 − t1) accounting for any dose or covariate events acting in the time interval. The equations were coded in the R language and used to simulate the time-course of drug amounts in each compartment of the systems. The equations were validated against commercial software [NONMEM (Beal, Sheiner, Boeckmann, & Bauer, 2009)] output to assess their capability to handle both complex dosage regimens and the effect of changes in covariate/ parameter values that may occur at discrete times within or between dosing intervals. Results: For all tested pharmacokinetic models, the time-course of drug amounts using the ADVAN-style analytical solutions were identical to NONMEM outputs to at least four significant figures, confirming the validity of the presented equations. Discussion: To our knowledge, this paper presents the ADVAN-style equations for common pharmacokinetic models in the literature for the first time. The presented ADVAN-style equations overcome obstacles to implementing the classical analytical solutions in software, and have speed advantages over solutions using differential equation solvers. The equations presented in this paper fill a gap in the pharmacokinetic literature, and it is expected that these equations will facilitate the investigation of useful open-source software for modelling pharmacokinetic data. © 2015 Elsevier Inc. All rights reserved.
1. Introduction The implementation of pharmacokinetic models in modelling software typically uses models expressed as either differential equations or as analytical solutions to compartmental systems. Analytical solutions are often many times faster than differential solutions, and are generally preferred for this reason when a choice is possible. The analytical solutions to 1, 2 and 3 compartment pharmacokinetic models have been known for many decades and are published in a number of sources (Rescigno & Segre, 1966; Wagner, 1975). However, when implementing these equations in programming languages such as R [www.r-project. org] for stochastic simulations or Bayes forecasting (Mould, Upton, & Wojciechowski, 2014) or in hierarchical modelling environments such as rccpbugs [https://github.com/armstrtw/rcppbugs] for population modelling, some limitations of the published forms of the equations are apparent. First, the representation of the effect of discrete changes in covariate/parameter values for different time intervals is not possible, ⁎ Corresponding author. Fax: +61 8 8302 2389. E-mail address: [email protected] (R.N. Upton).
http://dx.doi.org/10.1016/j.vascn.2015.03.004 1056-8719/© 2015 Elsevier Inc. All rights reserved.
and second the representation of complex dose regimens would require the principle of superposition (Thron, 1974), which is difficult and clumsy to code. Both of these limitations may be a significant obstacle to implementing open-source pharmacokinetic software that is broadly useful in its application. The proprietary population modelling software NONMEM (Beal, Sheiner, Boeckmann, & Bauer, 2009) implements analytical solutions to common compartment models via its ADVAN routines (e.g. ADVAN1 for a one compartment linear model, ADVAN3 for a two compartment linear model and ADVAN11 for a three compartment linear model). These routines are able to handle both complex dose regimens and the effect of discrete changes in covariate values that may occur at discrete times within or between dosing intervals by “advancing” the kinetic system from one state time point to the next (Beal et al., 2009). The use of this type of “ADVAN” solution requires recasting the analytical equations for common pharmacokinetic models in a form where they can be used to calculate the amount of drug in each compartment of the system at time t2 given the amount in each compartment at time t1 as a starting point, together with any dose or covariate factors acting in the period t2 − t1. This report presents the analytical solutions to the
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
43
ADVAN-style equations in the public domain for the first time. It is expected that this will facilitate the investigation of useful open-source software for modelling and simulating pharmacokinetic data. The aims of this paper are: 1 Solve and report the ADVAN-style equations for 1, 2 and 3 compartment pharmacokinetic linear models for intravenous (IV) and firstorder absorption drug administration. 2 Compare the output of the equations with the output of NONMEM for the same models. 3 Show example code for the implementation of such equations in the R language. 2. Methods 2.1. ADVAN-style equations Solutions of the ADVAN-style equations for 1, 2, and 3 compartment pharmacokinetic linear models were derived using Laplace transforms for finding analytical solutions of compartmental models (Rescigno & Segre, 1966; Wagner, 1975). All state variables were expressed as amounts (X) rather than concentrations to ensure mass balance. The solution procedure consisted of the following steps: (1) expressing the mass transfer of a drug for each pharmacokinetic model through a system of differential equations with the initial conditions included as part of the solution process. (2) Taking Laplace transforms of each term in the differential equations (i.e. transform equations from the t-domain into the s-domain). (3) Solving the Laplace of state variable(s) of interest. (4) Finally, taking the inverse of the Laplace of the state variable(s) (i.e. transform equations from the s-domain back into the t-domain). The solution procedure is depicted in Fig. 1. The Laplace transform of a function F(t) is defined by Z Lf F ðt Þg ¼ f ðsÞ ¼
∞
−st
e
F ðt Þdt:
ð1Þ
0
Fig. 2. Schematic representations of 1, 2, and 3-compartment pharmacokinetic models with IV bolus administration (Dose = D) expressed in terms of the amount of drug in each compartment. X1, X2, or X3 represents a drug amount in the corresponding compartment at the end of a time period (t) where (t = t2 − t1). A1, A2, or A3 represents the initial drug amount in the corresponding compartment at the beginning of a time interval (t). The transfer into and out of the compartments is governed by the micro-rate constants (time−1) given the symbol k (e.g. k12 between the first and second compartment, k13 between the first and third compartment).
2.1.1. Intravenous pharmacokinetic models 2.1.1.1. IV bolus administration. A one compartment IV bolus is the simplest pharmacokinetic model and is parameterized as presented in (Fig. 2a). The differential equation that describes the model is: dX 1 ¼ −k10 X 1 dt
ð2Þ
where X1 is a drug amount in the central compartment and k10 is the first-order elimination rate constant. Laplace transform of the variable of interest (X1) is taken over both sides of the equation: sLfX 1 g−X 1 ð0Þ ¼ −k10 LfX 1 g
Fig. 1. The steps involved in solving differential equations using Laplace transform. The method converts a complex differential equation into a simple algebraic equation that can be solved easily. The dashed line indicates that the differential equation may be solved by integration. However, integration might be difficult for complex models.
44
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
where X1(0) represents the initial drug amount in the compartment. In the case of ADVAN-style solutions, it represents the initial drug amount at the beginning of each time interval (t) where (t = t2 − t1) and is denoted as (A) in Fig. 2. Therefore, L{X1} can be written as: LfX 1 g ¼ X 1 ðsÞ ¼
A1 ðs þ k10 Þ
ð3Þ
−k10 t
X 1 ðt Þ ¼ A1 e
ð4Þ
where X1(t) is the drug amount at the end of a time interval (t), A1 is the initial drug amount in the compartment at the beginning of each time interval, and k10 is the first-order elimination rate constant. The differential equations describing a 2-compartment IV bolus model (Fig. 2b) are: dX 1 ¼ k21 X 2 − E1 X 1 dt dX 2 ¼ k12 X 1 − E2 X 2 dt
ð5Þ
where X1, and X2 are the drug amounts in the central and peripheral compartments, respectively. E1 and E2 are the sums of the elimination rate constants from the central, and peripheral compartments, respectively (E1 = k10 + k12 and E2 = k20 + k21). Taking Laplace transform of Eq. (5) and writing the coefficients in matrix notion gives: ðs þ E1 Þ −k12
−k21 ðs þ E2 Þ
A1 X 1 ðsÞ ¼ X 2 ðsÞ A2
ð6Þ
the amount of drug in the central and peripheral compartments in the Laplace domain (X1(s), and X2(s), respectively) is given by the ratio of the determinants of two matrices: A −k21 1 A ðs þ E Þ 2 2 X 1 ðsÞ ¼ ð −k21 s þ E1 Þ −k ðs þ E 2 Þ 12 ðs þ E Þ A 1 1 −k A2 12 X 2 ðsÞ ¼ −k21 ðs þ E1 Þ : −k ðs þ E Þ 12
ð7Þ
2
The solution to lower matrix evaluates to a quadratic polynomial in s: 2
s þ ðE1 þ E2 Þs þ E1 E2 −k12 k21
ð8Þ
which can be written for convenience as: ðs þ λ1 Þðs þ λ2 Þ
ð9Þ
where the negative values of λ1 and λ2 are the real roots of the quadratic polynomial of Eq. (8) and represent the hybrid-rate constants for the 2compartment model. Substituting Eq. (9) into Eq. (7) gives:
ð10Þ
The inverse Laplace transform of Eq. (10) gives: −λ1 t
½ðA1 E2 þ A2 k21 Þ−A1 λ1 e
ð12Þ
Applying a similar approach for the 3-compartment IV bolus model (Fig. 2c), the amount of drug in each compartment in the Laplace domain is given by the ratio of the determinants as described in Eq. (13). A −k12 −k31 1 A ðs þ E Þ 0 2 2 A 0 ð s þ E Þ 3 3 X 1 ðsÞ ¼ ð Þ −k −k s þ E 1 12 31 −k ðs þ E2 Þ 0 12 −k 0 ð s þ E Þ 31 3 ðs þ E Þ A −k 1 1 31 −k A2 0 12 −k A ð s þ E Þ 31 3 3 X 2 ðsÞ ¼ ð Þ −k −k s þ E 1 12 31 −k ðs þ E2 Þ 0 12 −k 0 ð s þ E Þ 31 3 ðs þ E Þ −k A 1 12 1 −k ðs þ E2 Þ A2 12 −k 0 A3 31 X 3 ðsÞ ¼ ð Þ −k −k s þ E 1 12 31 −k ðs þ E2 Þ 0 12 −k 0 ðs þ E3 Þ 31
ð13Þ
where E1, E2, and E3 are the sum of the elimination rate constants from the central, first-peripheral, and second-peripheral compartments, respectively (E1 = k12 + k13 + k10, E2 = k21 + k20, and E3 = k31 + k30). A1, A2 and A3 are the initial drug amounts in the central, firstperipheral, and second-peripheral compartments, respectively, at the beginning of each time interval (t). The solution to the lower matrix evaluates to a cubic polynomial in s: 3
2
s þ as þ bs þ c
ð14Þ
where: a ¼ E1 þ E2 þ E3 b ¼ E1 E2 þ E3 ðE1 þ E2 Þ−k12 k21 −k13 k31 c ¼ E1 E2 E3 −E3 k12 k21 −E2 k13 k31 :
ð15Þ
The cubic polynomial of Eq. (14) can be written for convenience as: ðs þ λ1 Þðs þ λ2 Þðs þ λ3 Þ
A ðs þ E2 Þ þ A2 k21 X 1 ðsÞ ¼ 1 ðs þ λ1 Þðs þ λ2 Þ A2 ðs þ E1 Þ þ A1 k12 : X 2 ðsÞ ¼ ðs þ λ1 Þðs þ λ2 Þ
ð16Þ
where the negative values of λ1, λ2, and λ3 are the real roots of the polynomial and the hybrid rate constants of the 3-compartment model. Substituting Eq. (16) in Eq. (13) gives: A1 ðs þ E2 Þðs þ E3 Þ þ A2 k21 ðs þ E3 Þ þ A3 k31 ðs þ E2 Þ ðs þ λ1 Þðs þ λ2 Þðs þ λ2 Þ A2 ½ðs þ E1 Þðs þ E3 Þ−k13 k31 þ A1 k12 ðs þ E3 Þ þ A3 k21 k31 X 2 ðsÞ ¼ ðs þ λ1 Þðs þ λ2 Þðs þ λ2 Þ A3 ½ðs þ E1 Þðs þ E2 Þ−k12 k21 þ A1 k13 ðs þ E2 Þ þ A2 k21 k13 X 3 ðsÞ ¼ ðs þ λ1 Þðs þ λ2 Þðs þ λ2 Þ
X 1 ðsÞ ¼ −λ t
−½ðA1 E2 þ A2 k21 Þ−A1 λ2 e 2 ðλ2 −λ1 Þ −λ t −λ t ð11Þ ½ðA2 E1 þ A1 k12 Þ−A2 λ1 e 1 −½ðA2 E1 þ A1 k12 Þ−A2 λ2 e 2 X 2 ðt Þ ¼ ðλ2 −λ1 Þ X 1 ðt Þ ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðE1 þ E2 Þ þ ðE1 þ E2 Þ2 −4ðE1 E2 −k12 k21 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 λ2 ¼ ðE1 þ E2 Þ− ðE1 þ E2 Þ2 −4ðE1 E2 −k12 k21 : 2 λ1 ¼
and the inverse Laplace transform is:
where X1(t), and X2(t) describe the time course of drug amounts in the central and peripheral compartments at the end of each time interval (t). A1, and A2 are the initial drug amounts at the beginning of each time interval. The hybrid rate constants (λ1, and λ2) are easily obtained from the micro-rate constants by Eq. (12).
ð17Þ
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
constants are calculated as described by Eq. (25). The same method is applied to obtain micro-rate constants for other pharmacokinetic models.
and the inverse Laplace of Eq. (17) are: X 1 ðt Þ ¼ A1
ðE2 −λ1 ÞðE3 −λ1 Þ −λ1 t ðE2 −λ2 ÞðE3 −λ2 Þ −λ2 t e e þ ðλ2 −λ1 Þðλ3 −λ1 Þ ðλ1 −λ2 Þðλ3 −λ2 Þ
ðE2 −λ3 ÞðE3 −λ3 Þ −λ3 t ðC−Bλ1 Þ −λ t e e 1 þ ðλ1 −λ3 Þðλ2 −λ3 Þ ðλ1 −λ2 Þðλ1 −λ3 Þ ðBλ2 −C Þ ðBλ3 −C Þ −λ t −λ t e 2 þ e 3 þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðλ1 −λ3 Þðλ3 −λ2 Þ
þ
X 2 ðt Þ ¼ A2
ðE1 −λ3 ÞðE3 −λ3 Þ −λ3 t ðI−A1 k12 λ1 Þ −λ t e e 1 þ ðλ1 −λ3 Þðλ2 −λ3 Þ ðλ1 −λ2 Þðλ1 −λ3 Þ
þ
ðA1 k12 λ2 −IÞ ðA1 k12 λ3 −IÞ −λ t −λ t e 2 þ e 3 ðλ1 −λ2 Þðλ2 −λ3 Þ ðλ1 −λ3 Þðλ3 −λ2 Þ
X 3 ðt Þ ¼ A3
ð18Þ k10 ¼
ðE1 −λ1 ÞðE3 −λ1 Þ −λ1 t ðE1 −λ2 ÞðE3 −λ2 Þ −λ2 t e e þ ðλ2 −λ1 Þðλ3 −λ1 Þ ðλ1 −λ2 Þðλ3 −λ2 Þ
þ
ð19Þ
ðE1 −λ1 ÞðE2 −λ1 Þ −λ1 t ðE1 −λ2 ÞðE2 −λ2 Þ −λ2 t e e þ ðλ2 −λ1 Þðλ3 −λ1 Þ ðλ1 −λ2 Þðλ3 −λ2 Þ
þ
ðE1 −λ3 ÞðE2 −λ3 Þ −λ3 t ð J−A1 k13 λ1 Þ −λ t e e 1 þ ðλ1 −λ3 Þðλ2 −λ3 Þ ðλ1 −λ2 Þðλ1 −λ3 Þ
þ
ðA1 k13 λ2 −J Þ ðA1 k13 λ3 −J Þ −λ t −λ t e 2 þ e 3 ðλ1 −λ2 Þðλ2 −λ3 Þ ðλ1 −λ3 Þðλ3 −λ2 Þ
ð20Þ
ð21Þ
A standard method for calculating the hybrid rate constants (λ1, λ2, and λ3) from the micro-rate constants has been reported previously by Upton (2004). By this method:
m¼ n¼
3 2a3 −9ab þ 27c 2
Q¼
ð22Þ
2.1.1.2. IV infusion administration. The parameterization of the IV infusion models is similar to the IV bolus models (Fig. 2) with the input being an infusion with an infusion rate (R) rather than a bolus. The ADVAN-style equations of IV infusion models are similar to the IV bolus equations with the addition of an extra term that relates to the drug input being an infusion. The equations for 1, 2, and 3compartment infusion models with an infusion rate (R) are presented in Eqs. (26)-(31). For a 1-compartment IV infusion model: R −k t −k t 1−e 10 þ A1 e 10 k10
X 2 ðt Þ ¼
pffiffiffiffiffiffiffiffiffi −Q ; Q b0 n β¼ − 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi γ¼ β2 þ α 2 θ ¼ atan2ðβ; α Þ: α¼
ð23Þ
The hybrid rate constants are calculated from the values of α, γ, and θ using the following equations: pffiffiffi 1 a θ θ þ γ =3 cos þ 3 sin 3 3 pffiffiffi 3 1 a θ θ − 3 sin λ2 ¼ þ γ =3 cos 3 3 3 1 a θ : λ3 ¼ − 2γ =3 cos 3 3 λ1 ¼
ð24Þ
The micro-rate constants (time−1) that govern the transfer into and out of the compartments can be obtained from the clearance and the volume of distribution values of the reference pharmacokinetic model. For a typical 3-compartment IV bolus model (Fig. 2c), the micro-rate
½ðA1 E2 þ R þ A2 k21 Þ−A1 λ1 e
−λ1 t
−λ t
−½ðA1 E2 þ R þ A2 k21 Þ−A1 λ2 e 2 ðλ2 −λ1 Þ " # ð27Þ −λ t −λ t 1 e 1 e 2 þRE2 þ − λ1 λ2 λ1 ðλ1 −λ2 Þ λ2 ðλ1 −λ2 Þ
3
The exponential terms of the cubic polynomial are calculated as follows:
ð26Þ
where R is the infusion rate, X1(t) is the amount of drug in the central compartment at the end of each time interval (t), A1 is the initial drug amount at the beginning of each time interval (t), and k10 is the firstorder elimination rate constant. For a 2-compartment IV infusion model:
27
n m þ : 4 27
ð25Þ
where CL is the clearance of the drug from the central compartment and Q2, Q3 are the inter-compartmental clearances between the central and first-peripheral, and second-peripheral compartments, respectively. V1, V2, V3 are the volumes of distribution in the central and the two peripheral compartments, respectively.
X 1 ðt Þ ¼
3b−a2
CL Q Q Q Q ; k12 ¼ 2 ; k21 ¼ 2 ; k13 ¼ 3 ; k31 ¼ 3 V1 V1 V2 V1 V3
X 1 ðt Þ ¼
where: B ¼ A2 k21 þ A3 k31 C ¼ A2 k21 E3 þ A3 k31 E2 I ¼ A1 k12 E3 −A2 k13 k31 þ A3 k31 k12 J ¼ A1 k13 E2 −A3 k12 k21 þ A2 k13 k21 :
45
½ðA2 E1 þ A1 k12 Þ−A2 λ1 e−λ1 t −½ðA2 E1 þ A1 k12 Þ−A2 λ2 e−λ2 t ðλ2 −λ1 Þ " # ð28Þ −λ t −λ t 1 e 1 e 2 − þ þRk12 λ1 λ2 λ1 ðλ1 −λ2 Þ λ2 ðλ1 −λ2 Þ
where X1(t), X2(t) are the drug amounts in the central and peripheral compartments, respectively, at the end of a time interval (t). A1, A2 are the initial drug amounts at the beginning of each time interval. E1, and E2 are the sums of the micro-rate constants from the central and peripheral compartments, respectively. λ1, and λ2 are the macro-rate constants calculated by Eq. (12). For a 3-compartment IV infusion model: E2 E 3 ðE2 −λ1 ÞðE3 −λ1 Þ −λ1 t ðE2 −λ2 ÞðE3 −λ2 Þ −λ2 t X 1 ðt Þ ¼ R − − e e λ1 λ2 λ3 λ1 ðλ2 −λ1 Þðλ3 −λ1 Þ λ2 ðλ1 −λ2 Þðλ3 −λ2 Þ ðE2 −λ3 ÞðE3 −λ3 Þ −λ3 t ðE2 −λ1 ÞðE3 −λ1 Þ −λ1 t − e e þ A1 ðλ2 −λ1 Þðλ3 −λ1 Þ λ3 ðλ1 −λ3 Þðλ2 −λ3 Þ ðE −λ2 ÞðE3 −λ2 Þ −λ2 t ðE2 −λ3 ÞðE3 −λ3 Þ −λ3 t ð29Þ þ 2 þ e e ðλ1 −λ2 Þðλ3 −λ2 Þ ðλ1 −λ3 Þðλ2 −λ3 Þ ðC−Bλ1 Þ ðBλ2 −C Þ −λ t −λ t e 1 þ e 2 þ ðλ1 −λ2 Þðλ1 −λ3 Þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðBλ3 −C Þ −λ t þ e 3 ðλ1 −λ3 Þðλ3 −λ2 Þ
46
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
E3 ðE3 −λ1 Þ ðE3 −λ2 Þ −λ t −λ t − e 1− e 2 λ1 λ2 λ3 λ1 ðλ2 −λ1 Þðλ3 −λ1 Þ λ2 ðλ1 −λ2 Þðλ3 −λ2 Þ ðE3 −λ3 Þ ðE1 −λ1 ÞðE3 −λ1 Þ −λ1 t −λ t − e 3 þ A2 e ðλ2 −λ1 Þðλ3 −λ1 Þ λ3 ðλ1 −λ3 Þðλ2 −λ3 Þ ðE −λ2 ÞðE3 −λ2 Þ −λ2 t ðE1 −λ3 ÞðE3 −λ3 Þ −λ3 t ð30Þ e e þ 1 þ ðλ1 −λ2 Þðλ3 −λ2 Þ ðλ1 −λ3 Þðλ2 −λ3 Þ ðI−A1 k12 λ1 Þ ðA1 k12 λ2 −IÞ −λ t −λ t þ e 1 þ e 2 ðλ1 −λ2 Þðλ1 −λ3 Þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðA1 k12 λ3 −I Þ −λ t e 3 þ ðλ1 −λ3 Þðλ3 −λ2 Þ
X 2 ðt Þ ¼ Rk12
ðE2 −λ1 Þ ðE2 −λ2 Þ E2 −λ t −λ t e 1− e 2 − λ1 λ2 λ3 λ1 ðλ2 −λ1 Þðλ3 −λ1 Þ λ2 ðλ1 −λ2 Þðλ3 −λ2 Þ ðE2 −λ3 Þ ðE1 −λ1 ÞðE2 −λ1 Þ −λ1 t −λ t − e 3 þ A3 e ðλ2 −λ1 Þðλ3 −λ1 Þ λ3 ðλ1 −λ3 Þðλ2 −λ3 Þ ðE −λ2 ÞðE2 −λ2 Þ −λ2 t ðE1 −λ3 ÞðE2 −λ3 Þ −λ3 t ð31Þ e e þ 1 þ ðλ1 −λ2 Þðλ3 −λ2 Þ ðλ1 −λ3 Þðλ2 −λ3 Þ ð J−A1 k13 λ1 Þ ðA1 k13 λ2 − J Þ −λ t −λ t þ e 1 þ e 2 ðλ1 −λ2 Þðλ1 −λ3 Þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðA1 k13 λ3 −J Þ −λ t þ e 3 ðλ1 −λ3 Þðλ3 −λ2 Þ
X 3 ðt Þ ¼ Rk13
where B, C, I, and J are given by Eq. (21). λ1, λ2, and λ3 are the real roots of the cubic polynomial of Eq. (14) calculated by Eq. (24).
2.1.2. First-order absorption models A schematic representation of the parameterization of 1, 2, and 3 compartment first-order absorption models is presented in Fig. 3. The ADVAN-style equations for the first-order absorption models are similar to the IV bolus equations with the addition of an extra term that relates to the drug input being an oral. The equations for 1, 2, and 3 compartment first-order absorption models are presented in Eqs. (32)–(39). For a 1-compartment first-order absorption model (Fig. 3a): −ka t
X 1 ðt Þ ¼ A1 e
ð32Þ
where X1(t) is the bioavailable drug amount in the absorption compartment at the end of a time interval (t), A1 is the initial bioavailable drug amount in the absorption compartment at the beginning of each time interval, and ka is the first-order absorption rate constant. Eq. (32) applies to all first-order absorption models (1, 2, and 3-compartment models).
X 2 ðt Þ ¼
A1 ka −k20 t −k t −k t e −e a þ A2 e 20 ðka −k20 Þ
ð33Þ
where X2(t) is the amount of drug in the central compartment at the end of a time interval (t), A1 and A2 are the initial drug amounts in the absorption and central compartments, respectively, at the beginning of each time interval, and k20 is the first-order elimination rate constant of the central compartment. For a 2-compartment first-order absorption model (Fig. 3b):
ðE3 −ka Þ ðE3 −λ1 Þ −k t −λ t e a þ e 1 ðλ1 −ka Þðλ2 −ka Þ ðλ2 −λ1 Þðka −λ1 Þ ðE3 −λ2 Þ −λ t e 2 þ ðλ1 −λ2 Þðka −λ2 Þ
X 2 ðt Þ ¼ A1 ka
þ
½ðA2 E3 þ A3 k32 Þ−A2 λ1 e−λ1 t −½ðA2 E3 þ A3 k32 Þ−A2 λ2 e−λ2 t ðλ2 −λ1 Þ ð34Þ "
X 3 ðt Þ ¼ A1 ka k23
þ
−k t
−λ t
−λ t
e a e 1 e 2 þ þ ðλ1 −ka Þðλ2 −ka Þ ðλ2 −λ1 Þðka −λ1 Þ ðλ1 −λ2 Þðka −λ2 Þ
½ðA3 E2 þ A2 k23 Þ−A3 λ1 e
−λ1 t
−½ðA3 E2 þ A2 k23 Þ−A3 λ2 e ðλ2 −λ1 Þ
#
−λ2 t
ð35Þ where E2 and E3 are the sum of the elimination rate constants from the central and peripheral compartments, respectively (E2 = k23 + k20 and E3 = k32 + k30). A2 and A3 are the initial drug amounts in the central and peripheral compartments, respectively, at the beginning of each time interval. λ1, and λ2 are the macro-rate constants calculated by Eq. (36). qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðE2 þ E3 Þ þ ðE2 þ E3 Þ2 −4ðE2 E3 −k23 k32 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðE2 þ E3 Þ− ðE2 þ E3 Þ2 −4ðE2 E3 −k23 k32 : λ2 ¼ 2
λ1 ¼
ð36Þ
For a 3-compartment first-order absorption model (Fig. 3c):
ðE3 −ka ÞðE4 −ka Þ ðE3 −λ1 ÞðE4 −λ1 Þ −k t −λ t e a þ e 1 ðλ1 −ka Þðλ2 −ka Þðλ3 −ka Þ ðka −λ1 Þðλ2 −λ1 Þðλ3 −λ1 Þ ðE3 −λ2 ÞðE4 −λ2 Þ ðE3 −λ3 ÞðE4 −λ3 Þ −λ t −λ t þ e 2 þ e 3 ðka −λ2 Þðλ1 −λ2 Þðλ3 −λ2 Þ ðka −λ3 Þðλ1 −λ3 Þðλ2 −λ3 Þ ðE3 −λ1 ÞðE4 −λ1 Þ −λ1 t ðE3 −λ2 ÞðE4 −λ2 Þ −λ2 t ð37Þ þ A2 þ e e ðλ2 −λ1 Þðλ3 −λ1 Þ ðλ1 −λ2 Þðλ3 −λ2 Þ ðE −λ3 ÞðE4 −λ3 Þ −λ3 t ðC−Bλ1 Þ −λ t þ þ 3 e e 1 ðλ1 −λ3 Þðλ2 −λ3 Þ ðλ1 −λ2 Þðλ1 −λ3 Þ ðBλ2 −C Þ ðBλ3 −C Þ −λ t −λ t þ e 2 þ e 3 ðλ1 −λ2 Þðλ2 −λ3 Þ ðλ1 −λ3 Þðλ3 −λ2 Þ
X 2 ðt Þ ¼ A1 ka
Fig. 3. Schematic representations of 1, 2, and 3-compartment pharmacokinetic models with oral administration and first-order absorption process expressed in terms of the amount of drug in each compartment. X1, X2, or X3 represents a drug amount in the corresponding compartment at the end of a time period (t) where (t = t2 − t1). A1, A2, or A3 represents the initial drug amount in the corresponding compartment at the beginning of a time interval (t). The transfer into and out of the compartments is governed by the micro-rate constants (time−1) given the symbol k (e.g. k12 between the first and second compartment, k13 between the first and third compartment).
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
X 3 ðt Þ ¼ A1 ka k23
47
ðE4 −ka Þ −k t e a ðλ1 −ka Þðλ2 −ka Þðλ3 −ka Þ
ðE4 −λ1 Þ ðE4 −λ3 Þ −λ t −λ t e 1 þ e 3 ðka −λ1 Þðλ2 −λ1 Þðλ3 −λ1 Þ ðka −λ3 Þðλ1 −λ3 Þðλ2 −λ3 Þ ðE2 −λ1 ÞðE4 −λ1 Þ −λ1 t ðE2 −λ2 ÞðE4 −λ2 Þ −λ2 t ð38Þ e e þ A3 þ ðλ2 −λ1 Þðλ3 −λ1 Þ ðλ1 −λ2 Þðλ3 −λ2 Þ ðE −λ3 ÞðE4 −λ3 Þ −λ3 t ðI−A2 k23 λ1 Þ −λ t e e 1 þ þ 2 ðλ1 −λ3 Þðλ2 −λ3 Þ ðλ1 −λ2 Þðλ1 −λ3 Þ ðA2 k23 λ2 −IÞ ðA2 k23 λ3 −IÞ −λ t −λ t e 2 þ e 3 þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðλ1 −λ3 Þðλ3 −λ2 Þ
þ
X 4 ðt Þ ¼ A1 ka k24
ðE3 −ka Þ −k t e a ðλ1 −ka Þðλ2 −ka Þðλ3 −ka Þ
ðE3 −λ1 Þ ðE3 −λ2 Þ −λ t −λ t e 1 þ e 2 ðka −λ1 Þðλ2 −λ1 Þðλ3 −λ1 Þ ðka −λ2 Þðλ1 −λ2 Þðλ3 −λ2 Þ ðE3 −λ3 Þ ðE2 −λ1 ÞðE3 −λ1 Þ −λ1 t −λ t e 3 þ A4 e þ ðka −λ3 Þðλ1 −λ3 Þðλ2 −λ3 Þ ðλ2 −λ1 Þðλ3 −λ1 Þ ð39Þ ðE −λ2 ÞðE3 −λ2 Þ −λ2 t ðE2 −λ3 ÞðE3 −λ3 Þ −λ3 t þ 2 þ e e ðλ1 −λ2 Þðλ3 −λ2 Þ ðλ1 −λ3 Þðλ2 −λ3 Þ ð J−A2 k24 λ1 Þ ðA2 k24 λ2 − J Þ −λ1 t −λ t þ þ e e 2 ðλ1 −λ2 Þðλ1 −λ3 Þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðA2 k24 λ3 −J Þ −λ t þ e 3 ðλ1 −λ3 Þðλ3 −λ2 Þ
þ
where E2, E3, and E4 are the sums of the elimination rate constants from the central, first-peripheral, and second-peripheral compartments, respectively (E2 = k 23 + k 24 + k 20 , E 3 = k 32 + k 30, and E4 = k42 + k40). A2, A3 and A4 are the initial drug amounts in the central, first-peripheral, and second-peripheral compartments, respectively, at the beginning of each time interval. λ1, λ2 and λ3 are the macro-rate constants calculated by Eq. (24). B, C, I, and J terms are calculated using Eq. (40). B ¼ C ¼ I ¼ J ¼
A3 k32 þ A4 k42 A3 k32 E4 þ A4 k42 E3 A2 k23 E4 −A3 k24 k42 þ A4 k42 k23 A2 k24 E3 −A4 k23 k32 þ A3 k24 k32 :
ð40Þ
2.2. Validation The reported ADVAN-style analytical solutions for 1, 2 and 3 compartment linear pharmacokinetic models of IV and first-order absorption drug administration were coded in R language (R Core Team, 2014) and used to simulate the time-course of drug amount in each compartment of a reference pharmacokinetic model. The “R codes” are programmed in a way to be capable of handling both variable dose regimens and the effect of discrete changes in covariate values that may occur at discrete times within or between dosing intervals. An exaggerated change in creatinine clearance was used as an example to assess equations' capability of accounting for the effect of the discrete change in the covariate value on the timecourse of drug amounts in the different compartments of a reference pharmacokinetic system. Example R codes of the different pharmacokinetic models (1, 2, 3, compartments IV bolus, infusion, firstorder oral administration) using the ADVAN-style equations are provided in the supplementary material. The examples were set for two subjects only for demonstration purposes; however, the “*.csv” files can be set to contain as many subjects as needed in NONMEM format. Variable doses and variable dosing regimens can be used for each subject. The simulated R output using the derived ADVANstyle equations was compared to simulations performed by NONMEM software (Beal et al., 2009) Version VII level 2 (ICON Development Solutions, Ellicott City, MD) using the same population
Fig. 4. Simulations of a three compartment first-order absorption model with discrete changes in creatinine clearance (CLCR) values. The simulated dose was a 100 mg given twice daily. CLCR was deliberately changed from 100 ml/min (Time b 15 h) to a 30 ml/min (Time N =15 h) to show the capability of the equations in accounting for the effect of CLCR change on central clearance. The dotted symbols represent R-simulations of the time-course of drug amounts in the absorption (blue), central (red), first-peripheral (green), and second-peripheral (pink) compartments. The continuous black line represents NONMEM simulations of the same compartments. Both R and NONMEM simulations were identical. The dashed red line represents the simulated drug amount in the central compartment without incorporating the covariate effect.
pharmacokinetic parameters. Creatinine clearance was added as a covariate on systemic clearance using a power model (Eq. (41)). CLi ¼ CLpop ðCLCR =100Þ
ð41Þ
where CLi, CLpop, and CLCR are the individual clearance, population clearance, and creatinine clearance, respectively. Creatinine clearance was deliberately changed from 100 ml/min (Time b 15 h) to a 30 ml/min (Time N = 15 h) to show the capability of the equations in accounting for the discrete changes in covariate values. The speed of simulating using the ADVAN-style analytical solutions was compared to using differential equations in the R software. The R package ‘deSolve’, a differential equation solver (Soetaert, Petzoldt, & Setzer, 2010), was utilized for this purpose. A two compartment multiple-infusion model was used as an example to determine the speed advantage obtained by using the analytical solutions. Simulations were performed for various numbers of subjects in the range of 10–1000 subjects. 3. Results The output of R simulations of the time-course of drug amount using ADVAN-style equations were identical to NONMEM output for all tested Table 1 Ordinary differential equations (ODE) versus analytical solutions: comparisona of computation time. Number of subjectsb 10 50 100 200 500 1000
Analytical solutions (secondsc)
ODE (secondsc)
0.34 1.57 3.18 7.14 20.3 50.47
3.23 4.79 12.35 15.00 33.69 84.57
a Based on a 2-compartment intravenous infusion model with three infusion doses administered over an interval of 100 h. Evaluations were done at 1 hour time intervals. b Simulations were performed using R (version 3.0.2) installed on a Dell® Power Edge R910 server with 4 × 10 core Xeon 2.26 GHz processors and 256 GB of RAM running Windows 8 Server SP1 64-bit. c Mean computation time of three runs.
48
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
intravenous and oral pharmacokinetic linear models to at least 4 significant figures. As an illustration, Fig. 4 shows a simulated time-course of drug amounts in the absorption, central, and the two peripheral compartments of a 3-compartment first-order absorption pharmacokinetic model. ADVAN-style equations were capable of handling multiple daily dosing and accounted for the effect of creatinine clearance change on the time-course of drug amount in the different compartments. A full set of graphs comparing R and NONMEM outputs of the different intravenous and first-order absorption models can be generated by running the R codes provided in the supplementary material. The simulations using the ADVAN-style analytical solutions were faster than the ones achieved using the equivalent differential equations in R. Table 1 shows the mean computation times to simulate three infusion doses for a range of a number of subjects. The simulations using the ADVAN-style analytical solutions were completed in 10–64% of the time compared to using the equivalent differential equations. Note that the simulation times were associated with some intrinsic variability presumably due to the dynamic allocation of system resources to the simulation task. 4. Discussion The main purpose of this paper was to derive ADVAN-style analytical solutions of common intravenous and first-order absorption pharmacokinetic models, implement them in an open-source software and make them available to the public domain. To our knowledge, this paper reports the ADVAN-style analytical solutions in the literature for the first time. The ADVAN-style analytical solutions overcome major limitations of the classical analytical solutions of common pharmacokinetic models in the literature. When applied in the pharmacokinetic modelling software, the latter solutions lack the capability of accounting for the effect of discrete changes in covariate/parameter values that may occur between or within dosing intervals and they are difficult to code for complex dose regimens. ADVAN-style analytical solutions were derived here for 1, 2, and 3 compartment linear pharmacokinetic models of IV bolus, infusion, and first-order drug absorption. An executive table summary of the analytical solutions for the different compartment models is provided in the supplementary material. In all derived equations, elimination was considered to be occurring from the central and peripheral compartments to allow for a general derivation of the equations and to consider the case when observed data are available for both the central and peripheral compartments. Readers should be aware that the pharmacokinetic system will be unidentifiable if only the central compartment was observed in which case the elimination from the peripheral compartment(s) should be set to zero. The presented equations are capable of handling arbitrary dosing regimens and they “advance” the solution of the model from one time point to the next, allowing for any dose or covariate factors to be accounted for. The ADVAN-style equations of the different pharmacokinetic models, with exaggerated changes in creatinine clearance values, were validated against NONMEM and both outputs were identical. The presented equations can be used for simulating the drug amount in each compartment of a reference pharmacokinetic model. One immediate application is for the simulation of population models coded in R. These models can be incorporated into user reactive web applications using the Shiny package of R (Wojciechowski, Hopkins, & Upton, 2015). As a real time-interface to model simulations, speed is important in these web applications when simulating larger populations (e.g. 1000 subjects). Our experience has shown that the ADVAN style analytical solutions to be noticeably faster than differential equations in this context. Ultimately, the ADVAN style equations could be incorporated into an open-source mixed-effect modelling framework to estimate population parameters. Speed is always desirable here, and capacity to handle complex dose regimens and covariate structures as provided by the ADVAN style equations is important for software that is useful for “real world”
problems. We have made preliminary investigations of the use of the rccpbugs package of R for this purpose. Note that the presented equations can be used to account for categorical covariates (e.g. Fed versus Fasted). In such case, an appropriate covariate model must be defined in the R code to describe covariate effects on the population pharmacokinetic parameters. The equations presented in the paper can also be further expanded to include, for example, metabolite concentrations or hypothetical effect site concentrations for PK-PD models. In order to showcase this applicability, examples showing the analytical solutions for a 3-compartment IV infusion, and a 3-compartment first-order absorption and elimination parent drug models, linked to a 1compartment first-order formation and elimination metabolite model were derived and are provided in the supplementary materials. Limitations of the presented work include: (i) handling doses administered using different routes of administration (e.g. IV bolus plus oral administration at the same time) cannot be done using the presented analytical solutions, where, in such scenarios, a new set of equations need to be derived to account for both drug inputs; (ii) solving for analytical solutions using Laplace transforms will increase in complexity as the number of states in the pharmacokinetic system increases; (iii) like all computer software, the ADVAN-style functions presented here may have problems when the values of the parameters or variables force a floating point overflow or underflow condition, or floating point rounding error. Users are advised to be vigilant for such scenarios in the software they choose to use for implementing the equations. The ADVAN-style analytical solutions presented herein fill a gap in the pharmacokinetic literature and may help facilitate the investigation of useful open-source software for modelling pharmacokinetic data. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.vascn.2015.03.004. Contributors Ahmad Y. Abuhelwa developed the equations and code. David J.R. Foster contributed to the validation of the equations. Richard N. Upton conceived the project and contributed to the development of the equations and code. All the authors helped in the drafting of the manuscript and have read and approved the final manuscript. Acknowledgement The authors acknowledge that the Australian Centre for Pharmacometrics is an initiative of the Australian Government as part of the National Collaborative Research Infrastructure Strategy. The authors acknowledge that Ahmad Y. Abuhelwa is a PhD student receiving an Endeavour Scholarship funded by the Department of Education and Training of the Australian Government (Scholarship ID no. 4088). References Beal, S., Sheiner, L. B., Boeckmann, A., & Bauer, R. J. (2009). NONMEM user's guides, part V, (1989–2009), Icon Development Solutions, Ellicott City, MD, USA. R Core Team (2014). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Mould, D. R., Upton, R. N., & Wojciechowski, J. (2014). Dashboard systems: Implementing pharmacometrics from bench to bedside. The AAPS Journal, 16(5), 925–937. Rescigno, A., & Segre, G. (1966). Drug and tracer kinetics. Blaisdell Pub. Co. Soetaert, K., Petzoldt, T., & Setzer, R. W. (2010). Solving differential equations in R: Package deSolve. Journal of Statistical Software, 33. Thron, C. (1974). Linearity and superposition in pharmacokinetics. Pharmacological Reviews, 26(1), 3–31. Upton, R. N. (2004). Calculating the hybrid (macro) rate constants of a threecompartment mamillary pharmacokinetic model from known micro-rate constants. Journal of Pharmacological and Toxicological Methods, 49(1), 65–68. Wagner, J. G. (1975). Fundamentals of clinical pharmacokinetics. Vol. 295, Illinois: Drug Intelligence Publications Hamilton. Wojciechowski, J., Hopkins, A., & Upton, R. (2015). Interactive pharmacometric applications using R and the shiny package. CPT: Pharmacometrics & Systems Pharmacology, 4(3).
Contents lists available at ScienceDirect
Journal of Pharmacological and Toxicological Methods journal homepage: www.elsevier.com/locate/jpharmtox
Original article
ADVAN-style analytical solutions for common pharmacokinetic models Ahmad Y. Abuhelwa, David J.R. Foster, Richard N. Upton ⁎ Australian Centre for Pharmacometrics and Sansom Institute, School of Pharmacy and Medical Sciences, University of South Australia, SA 5000, Australia
a r t i c l e
i n f o
Article history: Received 30 October 2014 Received in revised form 25 March 2015 Accepted 27 March 2015 Available online 2 April 2015 Keywords: ADVAN Analytical solutions Intravenous Methods Modelling NONMEM Oral Pharmacokinetics R-language Simulations
a b s t r a c t Introduction: The analytical solutions to compartmental pharmacokinetic models are well known, but have not been presented in a form that easily allows for complex dosing regimen and changes in covariate/parameter values that may occur at discrete times within and/or between dosing intervals. Methods: Laplace transforms were used to derive ADVAN-style analytical solutions for 1, 2, and 3 compartment pharmacokinetic linear models of intravenous and first-order absorption drug administration. The equations calculate the change in drug amounts in each compartment of the model over a time interval (t; t = t2 − t1) accounting for any dose or covariate events acting in the time interval. The equations were coded in the R language and used to simulate the time-course of drug amounts in each compartment of the systems. The equations were validated against commercial software [NONMEM (Beal, Sheiner, Boeckmann, & Bauer, 2009)] output to assess their capability to handle both complex dosage regimens and the effect of changes in covariate/ parameter values that may occur at discrete times within or between dosing intervals. Results: For all tested pharmacokinetic models, the time-course of drug amounts using the ADVAN-style analytical solutions were identical to NONMEM outputs to at least four significant figures, confirming the validity of the presented equations. Discussion: To our knowledge, this paper presents the ADVAN-style equations for common pharmacokinetic models in the literature for the first time. The presented ADVAN-style equations overcome obstacles to implementing the classical analytical solutions in software, and have speed advantages over solutions using differential equation solvers. The equations presented in this paper fill a gap in the pharmacokinetic literature, and it is expected that these equations will facilitate the investigation of useful open-source software for modelling pharmacokinetic data. © 2015 Elsevier Inc. All rights reserved.
1. Introduction The implementation of pharmacokinetic models in modelling software typically uses models expressed as either differential equations or as analytical solutions to compartmental systems. Analytical solutions are often many times faster than differential solutions, and are generally preferred for this reason when a choice is possible. The analytical solutions to 1, 2 and 3 compartment pharmacokinetic models have been known for many decades and are published in a number of sources (Rescigno & Segre, 1966; Wagner, 1975). However, when implementing these equations in programming languages such as R [www.r-project. org] for stochastic simulations or Bayes forecasting (Mould, Upton, & Wojciechowski, 2014) or in hierarchical modelling environments such as rccpbugs [https://github.com/armstrtw/rcppbugs] for population modelling, some limitations of the published forms of the equations are apparent. First, the representation of the effect of discrete changes in covariate/parameter values for different time intervals is not possible, ⁎ Corresponding author. Fax: +61 8 8302 2389. E-mail address: [email protected] (R.N. Upton).
http://dx.doi.org/10.1016/j.vascn.2015.03.004 1056-8719/© 2015 Elsevier Inc. All rights reserved.
and second the representation of complex dose regimens would require the principle of superposition (Thron, 1974), which is difficult and clumsy to code. Both of these limitations may be a significant obstacle to implementing open-source pharmacokinetic software that is broadly useful in its application. The proprietary population modelling software NONMEM (Beal, Sheiner, Boeckmann, & Bauer, 2009) implements analytical solutions to common compartment models via its ADVAN routines (e.g. ADVAN1 for a one compartment linear model, ADVAN3 for a two compartment linear model and ADVAN11 for a three compartment linear model). These routines are able to handle both complex dose regimens and the effect of discrete changes in covariate values that may occur at discrete times within or between dosing intervals by “advancing” the kinetic system from one state time point to the next (Beal et al., 2009). The use of this type of “ADVAN” solution requires recasting the analytical equations for common pharmacokinetic models in a form where they can be used to calculate the amount of drug in each compartment of the system at time t2 given the amount in each compartment at time t1 as a starting point, together with any dose or covariate factors acting in the period t2 − t1. This report presents the analytical solutions to the
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
43
ADVAN-style equations in the public domain for the first time. It is expected that this will facilitate the investigation of useful open-source software for modelling and simulating pharmacokinetic data. The aims of this paper are: 1 Solve and report the ADVAN-style equations for 1, 2 and 3 compartment pharmacokinetic linear models for intravenous (IV) and firstorder absorption drug administration. 2 Compare the output of the equations with the output of NONMEM for the same models. 3 Show example code for the implementation of such equations in the R language. 2. Methods 2.1. ADVAN-style equations Solutions of the ADVAN-style equations for 1, 2, and 3 compartment pharmacokinetic linear models were derived using Laplace transforms for finding analytical solutions of compartmental models (Rescigno & Segre, 1966; Wagner, 1975). All state variables were expressed as amounts (X) rather than concentrations to ensure mass balance. The solution procedure consisted of the following steps: (1) expressing the mass transfer of a drug for each pharmacokinetic model through a system of differential equations with the initial conditions included as part of the solution process. (2) Taking Laplace transforms of each term in the differential equations (i.e. transform equations from the t-domain into the s-domain). (3) Solving the Laplace of state variable(s) of interest. (4) Finally, taking the inverse of the Laplace of the state variable(s) (i.e. transform equations from the s-domain back into the t-domain). The solution procedure is depicted in Fig. 1. The Laplace transform of a function F(t) is defined by Z Lf F ðt Þg ¼ f ðsÞ ¼
∞
−st
e
F ðt Þdt:
ð1Þ
0
Fig. 2. Schematic representations of 1, 2, and 3-compartment pharmacokinetic models with IV bolus administration (Dose = D) expressed in terms of the amount of drug in each compartment. X1, X2, or X3 represents a drug amount in the corresponding compartment at the end of a time period (t) where (t = t2 − t1). A1, A2, or A3 represents the initial drug amount in the corresponding compartment at the beginning of a time interval (t). The transfer into and out of the compartments is governed by the micro-rate constants (time−1) given the symbol k (e.g. k12 between the first and second compartment, k13 between the first and third compartment).
2.1.1. Intravenous pharmacokinetic models 2.1.1.1. IV bolus administration. A one compartment IV bolus is the simplest pharmacokinetic model and is parameterized as presented in (Fig. 2a). The differential equation that describes the model is: dX 1 ¼ −k10 X 1 dt
ð2Þ
where X1 is a drug amount in the central compartment and k10 is the first-order elimination rate constant. Laplace transform of the variable of interest (X1) is taken over both sides of the equation: sLfX 1 g−X 1 ð0Þ ¼ −k10 LfX 1 g
Fig. 1. The steps involved in solving differential equations using Laplace transform. The method converts a complex differential equation into a simple algebraic equation that can be solved easily. The dashed line indicates that the differential equation may be solved by integration. However, integration might be difficult for complex models.
44
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
where X1(0) represents the initial drug amount in the compartment. In the case of ADVAN-style solutions, it represents the initial drug amount at the beginning of each time interval (t) where (t = t2 − t1) and is denoted as (A) in Fig. 2. Therefore, L{X1} can be written as: LfX 1 g ¼ X 1 ðsÞ ¼
A1 ðs þ k10 Þ
ð3Þ
−k10 t
X 1 ðt Þ ¼ A1 e
ð4Þ
where X1(t) is the drug amount at the end of a time interval (t), A1 is the initial drug amount in the compartment at the beginning of each time interval, and k10 is the first-order elimination rate constant. The differential equations describing a 2-compartment IV bolus model (Fig. 2b) are: dX 1 ¼ k21 X 2 − E1 X 1 dt dX 2 ¼ k12 X 1 − E2 X 2 dt
ð5Þ
where X1, and X2 are the drug amounts in the central and peripheral compartments, respectively. E1 and E2 are the sums of the elimination rate constants from the central, and peripheral compartments, respectively (E1 = k10 + k12 and E2 = k20 + k21). Taking Laplace transform of Eq. (5) and writing the coefficients in matrix notion gives: ðs þ E1 Þ −k12
−k21 ðs þ E2 Þ
A1 X 1 ðsÞ ¼ X 2 ðsÞ A2
ð6Þ
the amount of drug in the central and peripheral compartments in the Laplace domain (X1(s), and X2(s), respectively) is given by the ratio of the determinants of two matrices: A −k21 1 A ðs þ E Þ 2 2 X 1 ðsÞ ¼ ð −k21 s þ E1 Þ −k ðs þ E 2 Þ 12 ðs þ E Þ A 1 1 −k A2 12 X 2 ðsÞ ¼ −k21 ðs þ E1 Þ : −k ðs þ E Þ 12
ð7Þ
2
The solution to lower matrix evaluates to a quadratic polynomial in s: 2
s þ ðE1 þ E2 Þs þ E1 E2 −k12 k21
ð8Þ
which can be written for convenience as: ðs þ λ1 Þðs þ λ2 Þ
ð9Þ
where the negative values of λ1 and λ2 are the real roots of the quadratic polynomial of Eq. (8) and represent the hybrid-rate constants for the 2compartment model. Substituting Eq. (9) into Eq. (7) gives:
ð10Þ
The inverse Laplace transform of Eq. (10) gives: −λ1 t
½ðA1 E2 þ A2 k21 Þ−A1 λ1 e
ð12Þ
Applying a similar approach for the 3-compartment IV bolus model (Fig. 2c), the amount of drug in each compartment in the Laplace domain is given by the ratio of the determinants as described in Eq. (13). A −k12 −k31 1 A ðs þ E Þ 0 2 2 A 0 ð s þ E Þ 3 3 X 1 ðsÞ ¼ ð Þ −k −k s þ E 1 12 31 −k ðs þ E2 Þ 0 12 −k 0 ð s þ E Þ 31 3 ðs þ E Þ A −k 1 1 31 −k A2 0 12 −k A ð s þ E Þ 31 3 3 X 2 ðsÞ ¼ ð Þ −k −k s þ E 1 12 31 −k ðs þ E2 Þ 0 12 −k 0 ð s þ E Þ 31 3 ðs þ E Þ −k A 1 12 1 −k ðs þ E2 Þ A2 12 −k 0 A3 31 X 3 ðsÞ ¼ ð Þ −k −k s þ E 1 12 31 −k ðs þ E2 Þ 0 12 −k 0 ðs þ E3 Þ 31
ð13Þ
where E1, E2, and E3 are the sum of the elimination rate constants from the central, first-peripheral, and second-peripheral compartments, respectively (E1 = k12 + k13 + k10, E2 = k21 + k20, and E3 = k31 + k30). A1, A2 and A3 are the initial drug amounts in the central, firstperipheral, and second-peripheral compartments, respectively, at the beginning of each time interval (t). The solution to the lower matrix evaluates to a cubic polynomial in s: 3
2
s þ as þ bs þ c
ð14Þ
where: a ¼ E1 þ E2 þ E3 b ¼ E1 E2 þ E3 ðE1 þ E2 Þ−k12 k21 −k13 k31 c ¼ E1 E2 E3 −E3 k12 k21 −E2 k13 k31 :
ð15Þ
The cubic polynomial of Eq. (14) can be written for convenience as: ðs þ λ1 Þðs þ λ2 Þðs þ λ3 Þ
A ðs þ E2 Þ þ A2 k21 X 1 ðsÞ ¼ 1 ðs þ λ1 Þðs þ λ2 Þ A2 ðs þ E1 Þ þ A1 k12 : X 2 ðsÞ ¼ ðs þ λ1 Þðs þ λ2 Þ
ð16Þ
where the negative values of λ1, λ2, and λ3 are the real roots of the polynomial and the hybrid rate constants of the 3-compartment model. Substituting Eq. (16) in Eq. (13) gives: A1 ðs þ E2 Þðs þ E3 Þ þ A2 k21 ðs þ E3 Þ þ A3 k31 ðs þ E2 Þ ðs þ λ1 Þðs þ λ2 Þðs þ λ2 Þ A2 ½ðs þ E1 Þðs þ E3 Þ−k13 k31 þ A1 k12 ðs þ E3 Þ þ A3 k21 k31 X 2 ðsÞ ¼ ðs þ λ1 Þðs þ λ2 Þðs þ λ2 Þ A3 ½ðs þ E1 Þðs þ E2 Þ−k12 k21 þ A1 k13 ðs þ E2 Þ þ A2 k21 k13 X 3 ðsÞ ¼ ðs þ λ1 Þðs þ λ2 Þðs þ λ2 Þ
X 1 ðsÞ ¼ −λ t
−½ðA1 E2 þ A2 k21 Þ−A1 λ2 e 2 ðλ2 −λ1 Þ −λ t −λ t ð11Þ ½ðA2 E1 þ A1 k12 Þ−A2 λ1 e 1 −½ðA2 E1 þ A1 k12 Þ−A2 λ2 e 2 X 2 ðt Þ ¼ ðλ2 −λ1 Þ X 1 ðt Þ ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðE1 þ E2 Þ þ ðE1 þ E2 Þ2 −4ðE1 E2 −k12 k21 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 λ2 ¼ ðE1 þ E2 Þ− ðE1 þ E2 Þ2 −4ðE1 E2 −k12 k21 : 2 λ1 ¼
and the inverse Laplace transform is:
where X1(t), and X2(t) describe the time course of drug amounts in the central and peripheral compartments at the end of each time interval (t). A1, and A2 are the initial drug amounts at the beginning of each time interval. The hybrid rate constants (λ1, and λ2) are easily obtained from the micro-rate constants by Eq. (12).
ð17Þ
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
constants are calculated as described by Eq. (25). The same method is applied to obtain micro-rate constants for other pharmacokinetic models.
and the inverse Laplace of Eq. (17) are: X 1 ðt Þ ¼ A1
ðE2 −λ1 ÞðE3 −λ1 Þ −λ1 t ðE2 −λ2 ÞðE3 −λ2 Þ −λ2 t e e þ ðλ2 −λ1 Þðλ3 −λ1 Þ ðλ1 −λ2 Þðλ3 −λ2 Þ
ðE2 −λ3 ÞðE3 −λ3 Þ −λ3 t ðC−Bλ1 Þ −λ t e e 1 þ ðλ1 −λ3 Þðλ2 −λ3 Þ ðλ1 −λ2 Þðλ1 −λ3 Þ ðBλ2 −C Þ ðBλ3 −C Þ −λ t −λ t e 2 þ e 3 þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðλ1 −λ3 Þðλ3 −λ2 Þ
þ
X 2 ðt Þ ¼ A2
ðE1 −λ3 ÞðE3 −λ3 Þ −λ3 t ðI−A1 k12 λ1 Þ −λ t e e 1 þ ðλ1 −λ3 Þðλ2 −λ3 Þ ðλ1 −λ2 Þðλ1 −λ3 Þ
þ
ðA1 k12 λ2 −IÞ ðA1 k12 λ3 −IÞ −λ t −λ t e 2 þ e 3 ðλ1 −λ2 Þðλ2 −λ3 Þ ðλ1 −λ3 Þðλ3 −λ2 Þ
X 3 ðt Þ ¼ A3
ð18Þ k10 ¼
ðE1 −λ1 ÞðE3 −λ1 Þ −λ1 t ðE1 −λ2 ÞðE3 −λ2 Þ −λ2 t e e þ ðλ2 −λ1 Þðλ3 −λ1 Þ ðλ1 −λ2 Þðλ3 −λ2 Þ
þ
ð19Þ
ðE1 −λ1 ÞðE2 −λ1 Þ −λ1 t ðE1 −λ2 ÞðE2 −λ2 Þ −λ2 t e e þ ðλ2 −λ1 Þðλ3 −λ1 Þ ðλ1 −λ2 Þðλ3 −λ2 Þ
þ
ðE1 −λ3 ÞðE2 −λ3 Þ −λ3 t ð J−A1 k13 λ1 Þ −λ t e e 1 þ ðλ1 −λ3 Þðλ2 −λ3 Þ ðλ1 −λ2 Þðλ1 −λ3 Þ
þ
ðA1 k13 λ2 −J Þ ðA1 k13 λ3 −J Þ −λ t −λ t e 2 þ e 3 ðλ1 −λ2 Þðλ2 −λ3 Þ ðλ1 −λ3 Þðλ3 −λ2 Þ
ð20Þ
ð21Þ
A standard method for calculating the hybrid rate constants (λ1, λ2, and λ3) from the micro-rate constants has been reported previously by Upton (2004). By this method:
m¼ n¼
3 2a3 −9ab þ 27c 2
Q¼
ð22Þ
2.1.1.2. IV infusion administration. The parameterization of the IV infusion models is similar to the IV bolus models (Fig. 2) with the input being an infusion with an infusion rate (R) rather than a bolus. The ADVAN-style equations of IV infusion models are similar to the IV bolus equations with the addition of an extra term that relates to the drug input being an infusion. The equations for 1, 2, and 3compartment infusion models with an infusion rate (R) are presented in Eqs. (26)-(31). For a 1-compartment IV infusion model: R −k t −k t 1−e 10 þ A1 e 10 k10
X 2 ðt Þ ¼
pffiffiffiffiffiffiffiffiffi −Q ; Q b0 n β¼ − 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi γ¼ β2 þ α 2 θ ¼ atan2ðβ; α Þ: α¼
ð23Þ
The hybrid rate constants are calculated from the values of α, γ, and θ using the following equations: pffiffiffi 1 a θ θ þ γ =3 cos þ 3 sin 3 3 pffiffiffi 3 1 a θ θ − 3 sin λ2 ¼ þ γ =3 cos 3 3 3 1 a θ : λ3 ¼ − 2γ =3 cos 3 3 λ1 ¼
ð24Þ
The micro-rate constants (time−1) that govern the transfer into and out of the compartments can be obtained from the clearance and the volume of distribution values of the reference pharmacokinetic model. For a typical 3-compartment IV bolus model (Fig. 2c), the micro-rate
½ðA1 E2 þ R þ A2 k21 Þ−A1 λ1 e
−λ1 t
−λ t
−½ðA1 E2 þ R þ A2 k21 Þ−A1 λ2 e 2 ðλ2 −λ1 Þ " # ð27Þ −λ t −λ t 1 e 1 e 2 þRE2 þ − λ1 λ2 λ1 ðλ1 −λ2 Þ λ2 ðλ1 −λ2 Þ
3
The exponential terms of the cubic polynomial are calculated as follows:
ð26Þ
where R is the infusion rate, X1(t) is the amount of drug in the central compartment at the end of each time interval (t), A1 is the initial drug amount at the beginning of each time interval (t), and k10 is the firstorder elimination rate constant. For a 2-compartment IV infusion model:
27
n m þ : 4 27
ð25Þ
where CL is the clearance of the drug from the central compartment and Q2, Q3 are the inter-compartmental clearances between the central and first-peripheral, and second-peripheral compartments, respectively. V1, V2, V3 are the volumes of distribution in the central and the two peripheral compartments, respectively.
X 1 ðt Þ ¼
3b−a2
CL Q Q Q Q ; k12 ¼ 2 ; k21 ¼ 2 ; k13 ¼ 3 ; k31 ¼ 3 V1 V1 V2 V1 V3
X 1 ðt Þ ¼
where: B ¼ A2 k21 þ A3 k31 C ¼ A2 k21 E3 þ A3 k31 E2 I ¼ A1 k12 E3 −A2 k13 k31 þ A3 k31 k12 J ¼ A1 k13 E2 −A3 k12 k21 þ A2 k13 k21 :
45
½ðA2 E1 þ A1 k12 Þ−A2 λ1 e−λ1 t −½ðA2 E1 þ A1 k12 Þ−A2 λ2 e−λ2 t ðλ2 −λ1 Þ " # ð28Þ −λ t −λ t 1 e 1 e 2 − þ þRk12 λ1 λ2 λ1 ðλ1 −λ2 Þ λ2 ðλ1 −λ2 Þ
where X1(t), X2(t) are the drug amounts in the central and peripheral compartments, respectively, at the end of a time interval (t). A1, A2 are the initial drug amounts at the beginning of each time interval. E1, and E2 are the sums of the micro-rate constants from the central and peripheral compartments, respectively. λ1, and λ2 are the macro-rate constants calculated by Eq. (12). For a 3-compartment IV infusion model: E2 E 3 ðE2 −λ1 ÞðE3 −λ1 Þ −λ1 t ðE2 −λ2 ÞðE3 −λ2 Þ −λ2 t X 1 ðt Þ ¼ R − − e e λ1 λ2 λ3 λ1 ðλ2 −λ1 Þðλ3 −λ1 Þ λ2 ðλ1 −λ2 Þðλ3 −λ2 Þ ðE2 −λ3 ÞðE3 −λ3 Þ −λ3 t ðE2 −λ1 ÞðE3 −λ1 Þ −λ1 t − e e þ A1 ðλ2 −λ1 Þðλ3 −λ1 Þ λ3 ðλ1 −λ3 Þðλ2 −λ3 Þ ðE −λ2 ÞðE3 −λ2 Þ −λ2 t ðE2 −λ3 ÞðE3 −λ3 Þ −λ3 t ð29Þ þ 2 þ e e ðλ1 −λ2 Þðλ3 −λ2 Þ ðλ1 −λ3 Þðλ2 −λ3 Þ ðC−Bλ1 Þ ðBλ2 −C Þ −λ t −λ t e 1 þ e 2 þ ðλ1 −λ2 Þðλ1 −λ3 Þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðBλ3 −C Þ −λ t þ e 3 ðλ1 −λ3 Þðλ3 −λ2 Þ
46
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
E3 ðE3 −λ1 Þ ðE3 −λ2 Þ −λ t −λ t − e 1− e 2 λ1 λ2 λ3 λ1 ðλ2 −λ1 Þðλ3 −λ1 Þ λ2 ðλ1 −λ2 Þðλ3 −λ2 Þ ðE3 −λ3 Þ ðE1 −λ1 ÞðE3 −λ1 Þ −λ1 t −λ t − e 3 þ A2 e ðλ2 −λ1 Þðλ3 −λ1 Þ λ3 ðλ1 −λ3 Þðλ2 −λ3 Þ ðE −λ2 ÞðE3 −λ2 Þ −λ2 t ðE1 −λ3 ÞðE3 −λ3 Þ −λ3 t ð30Þ e e þ 1 þ ðλ1 −λ2 Þðλ3 −λ2 Þ ðλ1 −λ3 Þðλ2 −λ3 Þ ðI−A1 k12 λ1 Þ ðA1 k12 λ2 −IÞ −λ t −λ t þ e 1 þ e 2 ðλ1 −λ2 Þðλ1 −λ3 Þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðA1 k12 λ3 −I Þ −λ t e 3 þ ðλ1 −λ3 Þðλ3 −λ2 Þ
X 2 ðt Þ ¼ Rk12
ðE2 −λ1 Þ ðE2 −λ2 Þ E2 −λ t −λ t e 1− e 2 − λ1 λ2 λ3 λ1 ðλ2 −λ1 Þðλ3 −λ1 Þ λ2 ðλ1 −λ2 Þðλ3 −λ2 Þ ðE2 −λ3 Þ ðE1 −λ1 ÞðE2 −λ1 Þ −λ1 t −λ t − e 3 þ A3 e ðλ2 −λ1 Þðλ3 −λ1 Þ λ3 ðλ1 −λ3 Þðλ2 −λ3 Þ ðE −λ2 ÞðE2 −λ2 Þ −λ2 t ðE1 −λ3 ÞðE2 −λ3 Þ −λ3 t ð31Þ e e þ 1 þ ðλ1 −λ2 Þðλ3 −λ2 Þ ðλ1 −λ3 Þðλ2 −λ3 Þ ð J−A1 k13 λ1 Þ ðA1 k13 λ2 − J Þ −λ t −λ t þ e 1 þ e 2 ðλ1 −λ2 Þðλ1 −λ3 Þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðA1 k13 λ3 −J Þ −λ t þ e 3 ðλ1 −λ3 Þðλ3 −λ2 Þ
X 3 ðt Þ ¼ Rk13
where B, C, I, and J are given by Eq. (21). λ1, λ2, and λ3 are the real roots of the cubic polynomial of Eq. (14) calculated by Eq. (24).
2.1.2. First-order absorption models A schematic representation of the parameterization of 1, 2, and 3 compartment first-order absorption models is presented in Fig. 3. The ADVAN-style equations for the first-order absorption models are similar to the IV bolus equations with the addition of an extra term that relates to the drug input being an oral. The equations for 1, 2, and 3 compartment first-order absorption models are presented in Eqs. (32)–(39). For a 1-compartment first-order absorption model (Fig. 3a): −ka t
X 1 ðt Þ ¼ A1 e
ð32Þ
where X1(t) is the bioavailable drug amount in the absorption compartment at the end of a time interval (t), A1 is the initial bioavailable drug amount in the absorption compartment at the beginning of each time interval, and ka is the first-order absorption rate constant. Eq. (32) applies to all first-order absorption models (1, 2, and 3-compartment models).
X 2 ðt Þ ¼
A1 ka −k20 t −k t −k t e −e a þ A2 e 20 ðka −k20 Þ
ð33Þ
where X2(t) is the amount of drug in the central compartment at the end of a time interval (t), A1 and A2 are the initial drug amounts in the absorption and central compartments, respectively, at the beginning of each time interval, and k20 is the first-order elimination rate constant of the central compartment. For a 2-compartment first-order absorption model (Fig. 3b):
ðE3 −ka Þ ðE3 −λ1 Þ −k t −λ t e a þ e 1 ðλ1 −ka Þðλ2 −ka Þ ðλ2 −λ1 Þðka −λ1 Þ ðE3 −λ2 Þ −λ t e 2 þ ðλ1 −λ2 Þðka −λ2 Þ
X 2 ðt Þ ¼ A1 ka
þ
½ðA2 E3 þ A3 k32 Þ−A2 λ1 e−λ1 t −½ðA2 E3 þ A3 k32 Þ−A2 λ2 e−λ2 t ðλ2 −λ1 Þ ð34Þ "
X 3 ðt Þ ¼ A1 ka k23
þ
−k t
−λ t
−λ t
e a e 1 e 2 þ þ ðλ1 −ka Þðλ2 −ka Þ ðλ2 −λ1 Þðka −λ1 Þ ðλ1 −λ2 Þðka −λ2 Þ
½ðA3 E2 þ A2 k23 Þ−A3 λ1 e
−λ1 t
−½ðA3 E2 þ A2 k23 Þ−A3 λ2 e ðλ2 −λ1 Þ
#
−λ2 t
ð35Þ where E2 and E3 are the sum of the elimination rate constants from the central and peripheral compartments, respectively (E2 = k23 + k20 and E3 = k32 + k30). A2 and A3 are the initial drug amounts in the central and peripheral compartments, respectively, at the beginning of each time interval. λ1, and λ2 are the macro-rate constants calculated by Eq. (36). qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðE2 þ E3 Þ þ ðE2 þ E3 Þ2 −4ðE2 E3 −k23 k32 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðE2 þ E3 Þ− ðE2 þ E3 Þ2 −4ðE2 E3 −k23 k32 : λ2 ¼ 2
λ1 ¼
ð36Þ
For a 3-compartment first-order absorption model (Fig. 3c):
ðE3 −ka ÞðE4 −ka Þ ðE3 −λ1 ÞðE4 −λ1 Þ −k t −λ t e a þ e 1 ðλ1 −ka Þðλ2 −ka Þðλ3 −ka Þ ðka −λ1 Þðλ2 −λ1 Þðλ3 −λ1 Þ ðE3 −λ2 ÞðE4 −λ2 Þ ðE3 −λ3 ÞðE4 −λ3 Þ −λ t −λ t þ e 2 þ e 3 ðka −λ2 Þðλ1 −λ2 Þðλ3 −λ2 Þ ðka −λ3 Þðλ1 −λ3 Þðλ2 −λ3 Þ ðE3 −λ1 ÞðE4 −λ1 Þ −λ1 t ðE3 −λ2 ÞðE4 −λ2 Þ −λ2 t ð37Þ þ A2 þ e e ðλ2 −λ1 Þðλ3 −λ1 Þ ðλ1 −λ2 Þðλ3 −λ2 Þ ðE −λ3 ÞðE4 −λ3 Þ −λ3 t ðC−Bλ1 Þ −λ t þ þ 3 e e 1 ðλ1 −λ3 Þðλ2 −λ3 Þ ðλ1 −λ2 Þðλ1 −λ3 Þ ðBλ2 −C Þ ðBλ3 −C Þ −λ t −λ t þ e 2 þ e 3 ðλ1 −λ2 Þðλ2 −λ3 Þ ðλ1 −λ3 Þðλ3 −λ2 Þ
X 2 ðt Þ ¼ A1 ka
Fig. 3. Schematic representations of 1, 2, and 3-compartment pharmacokinetic models with oral administration and first-order absorption process expressed in terms of the amount of drug in each compartment. X1, X2, or X3 represents a drug amount in the corresponding compartment at the end of a time period (t) where (t = t2 − t1). A1, A2, or A3 represents the initial drug amount in the corresponding compartment at the beginning of a time interval (t). The transfer into and out of the compartments is governed by the micro-rate constants (time−1) given the symbol k (e.g. k12 between the first and second compartment, k13 between the first and third compartment).
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
X 3 ðt Þ ¼ A1 ka k23
47
ðE4 −ka Þ −k t e a ðλ1 −ka Þðλ2 −ka Þðλ3 −ka Þ
ðE4 −λ1 Þ ðE4 −λ3 Þ −λ t −λ t e 1 þ e 3 ðka −λ1 Þðλ2 −λ1 Þðλ3 −λ1 Þ ðka −λ3 Þðλ1 −λ3 Þðλ2 −λ3 Þ ðE2 −λ1 ÞðE4 −λ1 Þ −λ1 t ðE2 −λ2 ÞðE4 −λ2 Þ −λ2 t ð38Þ e e þ A3 þ ðλ2 −λ1 Þðλ3 −λ1 Þ ðλ1 −λ2 Þðλ3 −λ2 Þ ðE −λ3 ÞðE4 −λ3 Þ −λ3 t ðI−A2 k23 λ1 Þ −λ t e e 1 þ þ 2 ðλ1 −λ3 Þðλ2 −λ3 Þ ðλ1 −λ2 Þðλ1 −λ3 Þ ðA2 k23 λ2 −IÞ ðA2 k23 λ3 −IÞ −λ t −λ t e 2 þ e 3 þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðλ1 −λ3 Þðλ3 −λ2 Þ
þ
X 4 ðt Þ ¼ A1 ka k24
ðE3 −ka Þ −k t e a ðλ1 −ka Þðλ2 −ka Þðλ3 −ka Þ
ðE3 −λ1 Þ ðE3 −λ2 Þ −λ t −λ t e 1 þ e 2 ðka −λ1 Þðλ2 −λ1 Þðλ3 −λ1 Þ ðka −λ2 Þðλ1 −λ2 Þðλ3 −λ2 Þ ðE3 −λ3 Þ ðE2 −λ1 ÞðE3 −λ1 Þ −λ1 t −λ t e 3 þ A4 e þ ðka −λ3 Þðλ1 −λ3 Þðλ2 −λ3 Þ ðλ2 −λ1 Þðλ3 −λ1 Þ ð39Þ ðE −λ2 ÞðE3 −λ2 Þ −λ2 t ðE2 −λ3 ÞðE3 −λ3 Þ −λ3 t þ 2 þ e e ðλ1 −λ2 Þðλ3 −λ2 Þ ðλ1 −λ3 Þðλ2 −λ3 Þ ð J−A2 k24 λ1 Þ ðA2 k24 λ2 − J Þ −λ1 t −λ t þ þ e e 2 ðλ1 −λ2 Þðλ1 −λ3 Þ ðλ1 −λ2 Þðλ2 −λ3 Þ ðA2 k24 λ3 −J Þ −λ t þ e 3 ðλ1 −λ3 Þðλ3 −λ2 Þ
þ
where E2, E3, and E4 are the sums of the elimination rate constants from the central, first-peripheral, and second-peripheral compartments, respectively (E2 = k 23 + k 24 + k 20 , E 3 = k 32 + k 30, and E4 = k42 + k40). A2, A3 and A4 are the initial drug amounts in the central, first-peripheral, and second-peripheral compartments, respectively, at the beginning of each time interval. λ1, λ2 and λ3 are the macro-rate constants calculated by Eq. (24). B, C, I, and J terms are calculated using Eq. (40). B ¼ C ¼ I ¼ J ¼
A3 k32 þ A4 k42 A3 k32 E4 þ A4 k42 E3 A2 k23 E4 −A3 k24 k42 þ A4 k42 k23 A2 k24 E3 −A4 k23 k32 þ A3 k24 k32 :
ð40Þ
2.2. Validation The reported ADVAN-style analytical solutions for 1, 2 and 3 compartment linear pharmacokinetic models of IV and first-order absorption drug administration were coded in R language (R Core Team, 2014) and used to simulate the time-course of drug amount in each compartment of a reference pharmacokinetic model. The “R codes” are programmed in a way to be capable of handling both variable dose regimens and the effect of discrete changes in covariate values that may occur at discrete times within or between dosing intervals. An exaggerated change in creatinine clearance was used as an example to assess equations' capability of accounting for the effect of the discrete change in the covariate value on the timecourse of drug amounts in the different compartments of a reference pharmacokinetic system. Example R codes of the different pharmacokinetic models (1, 2, 3, compartments IV bolus, infusion, firstorder oral administration) using the ADVAN-style equations are provided in the supplementary material. The examples were set for two subjects only for demonstration purposes; however, the “*.csv” files can be set to contain as many subjects as needed in NONMEM format. Variable doses and variable dosing regimens can be used for each subject. The simulated R output using the derived ADVANstyle equations was compared to simulations performed by NONMEM software (Beal et al., 2009) Version VII level 2 (ICON Development Solutions, Ellicott City, MD) using the same population
Fig. 4. Simulations of a three compartment first-order absorption model with discrete changes in creatinine clearance (CLCR) values. The simulated dose was a 100 mg given twice daily. CLCR was deliberately changed from 100 ml/min (Time b 15 h) to a 30 ml/min (Time N =15 h) to show the capability of the equations in accounting for the effect of CLCR change on central clearance. The dotted symbols represent R-simulations of the time-course of drug amounts in the absorption (blue), central (red), first-peripheral (green), and second-peripheral (pink) compartments. The continuous black line represents NONMEM simulations of the same compartments. Both R and NONMEM simulations were identical. The dashed red line represents the simulated drug amount in the central compartment without incorporating the covariate effect.
pharmacokinetic parameters. Creatinine clearance was added as a covariate on systemic clearance using a power model (Eq. (41)). CLi ¼ CLpop ðCLCR =100Þ
ð41Þ
where CLi, CLpop, and CLCR are the individual clearance, population clearance, and creatinine clearance, respectively. Creatinine clearance was deliberately changed from 100 ml/min (Time b 15 h) to a 30 ml/min (Time N = 15 h) to show the capability of the equations in accounting for the discrete changes in covariate values. The speed of simulating using the ADVAN-style analytical solutions was compared to using differential equations in the R software. The R package ‘deSolve’, a differential equation solver (Soetaert, Petzoldt, & Setzer, 2010), was utilized for this purpose. A two compartment multiple-infusion model was used as an example to determine the speed advantage obtained by using the analytical solutions. Simulations were performed for various numbers of subjects in the range of 10–1000 subjects. 3. Results The output of R simulations of the time-course of drug amount using ADVAN-style equations were identical to NONMEM output for all tested Table 1 Ordinary differential equations (ODE) versus analytical solutions: comparisona of computation time. Number of subjectsb 10 50 100 200 500 1000
Analytical solutions (secondsc)
ODE (secondsc)
0.34 1.57 3.18 7.14 20.3 50.47
3.23 4.79 12.35 15.00 33.69 84.57
a Based on a 2-compartment intravenous infusion model with three infusion doses administered over an interval of 100 h. Evaluations were done at 1 hour time intervals. b Simulations were performed using R (version 3.0.2) installed on a Dell® Power Edge R910 server with 4 × 10 core Xeon 2.26 GHz processors and 256 GB of RAM running Windows 8 Server SP1 64-bit. c Mean computation time of three runs.
48
A.Y. Abuhelwa et al. / Journal of Pharmacological and Toxicological Methods 73 (2015) 42–48
intravenous and oral pharmacokinetic linear models to at least 4 significant figures. As an illustration, Fig. 4 shows a simulated time-course of drug amounts in the absorption, central, and the two peripheral compartments of a 3-compartment first-order absorption pharmacokinetic model. ADVAN-style equations were capable of handling multiple daily dosing and accounted for the effect of creatinine clearance change on the time-course of drug amount in the different compartments. A full set of graphs comparing R and NONMEM outputs of the different intravenous and first-order absorption models can be generated by running the R codes provided in the supplementary material. The simulations using the ADVAN-style analytical solutions were faster than the ones achieved using the equivalent differential equations in R. Table 1 shows the mean computation times to simulate three infusion doses for a range of a number of subjects. The simulations using the ADVAN-style analytical solutions were completed in 10–64% of the time compared to using the equivalent differential equations. Note that the simulation times were associated with some intrinsic variability presumably due to the dynamic allocation of system resources to the simulation task. 4. Discussion The main purpose of this paper was to derive ADVAN-style analytical solutions of common intravenous and first-order absorption pharmacokinetic models, implement them in an open-source software and make them available to the public domain. To our knowledge, this paper reports the ADVAN-style analytical solutions in the literature for the first time. The ADVAN-style analytical solutions overcome major limitations of the classical analytical solutions of common pharmacokinetic models in the literature. When applied in the pharmacokinetic modelling software, the latter solutions lack the capability of accounting for the effect of discrete changes in covariate/parameter values that may occur between or within dosing intervals and they are difficult to code for complex dose regimens. ADVAN-style analytical solutions were derived here for 1, 2, and 3 compartment linear pharmacokinetic models of IV bolus, infusion, and first-order drug absorption. An executive table summary of the analytical solutions for the different compartment models is provided in the supplementary material. In all derived equations, elimination was considered to be occurring from the central and peripheral compartments to allow for a general derivation of the equations and to consider the case when observed data are available for both the central and peripheral compartments. Readers should be aware that the pharmacokinetic system will be unidentifiable if only the central compartment was observed in which case the elimination from the peripheral compartment(s) should be set to zero. The presented equations are capable of handling arbitrary dosing regimens and they “advance” the solution of the model from one time point to the next, allowing for any dose or covariate factors to be accounted for. The ADVAN-style equations of the different pharmacokinetic models, with exaggerated changes in creatinine clearance values, were validated against NONMEM and both outputs were identical. The presented equations can be used for simulating the drug amount in each compartment of a reference pharmacokinetic model. One immediate application is for the simulation of population models coded in R. These models can be incorporated into user reactive web applications using the Shiny package of R (Wojciechowski, Hopkins, & Upton, 2015). As a real time-interface to model simulations, speed is important in these web applications when simulating larger populations (e.g. 1000 subjects). Our experience has shown that the ADVAN style analytical solutions to be noticeably faster than differential equations in this context. Ultimately, the ADVAN style equations could be incorporated into an open-source mixed-effect modelling framework to estimate population parameters. Speed is always desirable here, and capacity to handle complex dose regimens and covariate structures as provided by the ADVAN style equations is important for software that is useful for “real world”
problems. We have made preliminary investigations of the use of the rccpbugs package of R for this purpose. Note that the presented equations can be used to account for categorical covariates (e.g. Fed versus Fasted). In such case, an appropriate covariate model must be defined in the R code to describe covariate effects on the population pharmacokinetic parameters. The equations presented in the paper can also be further expanded to include, for example, metabolite concentrations or hypothetical effect site concentrations for PK-PD models. In order to showcase this applicability, examples showing the analytical solutions for a 3-compartment IV infusion, and a 3-compartment first-order absorption and elimination parent drug models, linked to a 1compartment first-order formation and elimination metabolite model were derived and are provided in the supplementary materials. Limitations of the presented work include: (i) handling doses administered using different routes of administration (e.g. IV bolus plus oral administration at the same time) cannot be done using the presented analytical solutions, where, in such scenarios, a new set of equations need to be derived to account for both drug inputs; (ii) solving for analytical solutions using Laplace transforms will increase in complexity as the number of states in the pharmacokinetic system increases; (iii) like all computer software, the ADVAN-style functions presented here may have problems when the values of the parameters or variables force a floating point overflow or underflow condition, or floating point rounding error. Users are advised to be vigilant for such scenarios in the software they choose to use for implementing the equations. The ADVAN-style analytical solutions presented herein fill a gap in the pharmacokinetic literature and may help facilitate the investigation of useful open-source software for modelling pharmacokinetic data. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.vascn.2015.03.004. Contributors Ahmad Y. Abuhelwa developed the equations and code. David J.R. Foster contributed to the validation of the equations. Richard N. Upton conceived the project and contributed to the development of the equations and code. All the authors helped in the drafting of the manuscript and have read and approved the final manuscript. Acknowledgement The authors acknowledge that the Australian Centre for Pharmacometrics is an initiative of the Australian Government as part of the National Collaborative Research Infrastructure Strategy. The authors acknowledge that Ahmad Y. Abuhelwa is a PhD student receiving an Endeavour Scholarship funded by the Department of Education and Training of the Australian Government (Scholarship ID no. 4088). References Beal, S., Sheiner, L. B., Boeckmann, A., & Bauer, R. J. (2009). NONMEM user's guides, part V, (1989–2009), Icon Development Solutions, Ellicott City, MD, USA. R Core Team (2014). R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Mould, D. R., Upton, R. N., & Wojciechowski, J. (2014). Dashboard systems: Implementing pharmacometrics from bench to bedside. The AAPS Journal, 16(5), 925–937. Rescigno, A., & Segre, G. (1966). Drug and tracer kinetics. Blaisdell Pub. Co. Soetaert, K., Petzoldt, T., & Setzer, R. W. (2010). Solving differential equations in R: Package deSolve. Journal of Statistical Software, 33. Thron, C. (1974). Linearity and superposition in pharmacokinetics. Pharmacological Reviews, 26(1), 3–31. Upton, R. N. (2004). Calculating the hybrid (macro) rate constants of a threecompartment mamillary pharmacokinetic model from known micro-rate constants. Journal of Pharmacological and Toxicological Methods, 49(1), 65–68. Wagner, J. G. (1975). Fundamentals of clinical pharmacokinetics. Vol. 295, Illinois: Drug Intelligence Publications Hamilton. Wojciechowski, J., Hopkins, A., & Upton, R. (2015). Interactive pharmacometric applications using R and the shiny package. CPT: Pharmacometrics & Systems Pharmacology, 4(3).
