J Pharmacokinet Pharmacodyn (2013) 40:639–650 DOI 10.1007/s10928-013-9336-y
ORIGINAL PAPER
Solution and implementation of distributed lifespan models Gilbert Koch • Johannes Schropp
Received: 25 July 2013 / Accepted: 8 October 2013 / Published online: 1 November 2013 Springer Science+Business Media New York 2013
Abstract We consider a population where every individual has a unique lifespan. After expiring of its lifespan the individual has to leave the population. A realistic approach to describe these lifespans is by a continuous distribution. Such distributed lifespan models (DLSMs) were introduced earlier in the indirect response context and consist of the mathematical convolution operator to describe the rate of change. Therefore, DLSMs could not directly be implemented in standard PKPD software. In this work we present the solution representation of DLSMs with and without a precursor population and an implementation strategy for DLSMs in ADAPT, NONMEM and MATLAB. We fit hemoglobin measurements from literature and investigate computational properties. Keywords Lifespan Distribution Convolution Population
After expiration of its lifespan the individual has to leave the population and maybe maturate into another population. A realistic approach to describe lifespans is by a continuous distribution. Such models were introduced by Krzyzanski et al. [1] to PKPD modeling in the indirect response context in 2006. In [1, 2] these models where applied to describe the stimulation of reticulocytes. In [3, 4] simpler models where every individual has the same lifespan are used to describe erythropoiesis. Let T be a random variable describing the distribution of the lifespan s of individuals based on a probability density function (PDF) l(s). The only assumption is that negative lifespans are not possible, more precisely l(s) = 0 for s \ 0. The expectation of the lifespan distribution is denoted by Z1 T ¼ E½T ¼ slðsÞds; 0
Introduction
the variance is
We consider the development of a population N(t) in time t where every individual has its own and unique lifespan s.
r2 ¼ Var½T ¼
Electronic supplementary material The online version of this article (doi:10.1007/s10928-013-9336-y) contains supplementary material, which is available to authorized users. G. Koch (&) Department of Pharmaceutical Sciences, State University of New York at Buffalo, 403 Kapoor Hall, Buffalo, NY 14214, USA e-mail: [email protected] J. Schropp Department of Mathematics and Statistics, Universita¨t Konstanz, PO Box 195, 78457 Konstanz, Germany
Z1
s2 lðsÞds T 2 ;
0
and the cumulative distribution function (CDF) is defined by LðsÞ ¼
Zs lðsÞds for
s 0:
0
In the special case of a point distribution (i.e. every individual has the same lifespan T) we have l(s) = d(s - T) where d is the Dirac delta function. A model describing the rate of change of a population reads
123
640
J Pharmacokinet Pharmacodyn (2013) 40:639–650
d NðtÞ ¼ kin ðtÞ kout ðtÞ; dt
Nð0Þ ¼ N 0 ;
ð1Þ
where the production (birth) rate is described by kin(t) and kout(t) denotes the loss (death) rate. The initial amount of individuals in the population at time point t = 0 is denoted by N0. Please note that only zero order terms appear in (1). The lifespan distribution is included into (1) via the convolution operator Z1 kout ðtÞ ¼ ðkin lÞðtÞ ¼ kin ðt sÞlðsÞds; ð2Þ
Fig. 1 The upper scheme is the basic distributed lifespan model, DLSM (5), (6) and the lower scheme shows the distributed lifespan model with precursor, DPLSM (7), (8)
0
see [1] for derivation. Hence, to calculate the outflow kout(t) at a time point t an integral over the lifespan interval I ¼ ½0; 1 has to be solved. Note that in (2) the past t \ 0 of the production kin(t) plays an important role and has to be defined. In this work we admit a non-constant past for the production kin(t) for t \ 0 which is an extension of traditional PKPD modeling and covers more pharmacological situations. Please note that (2) reduces to kout(t) = kin(t - T) in case of a point distributed lifespan, see [1]. The initial value N(0) in (1) describes the amount of born and already died individuals in the past up to the starting time point t = 0 and is calculated by
from that population maturate into a follow-up population M(t). Mathematically, this means that kin and lN from (5), (6) have to be replaced by kin * lN and lM. Hence, for a DPLSM we additionally obtain Z1 d MðtÞ ¼ ðkin lN ÞðtÞ ðkin lN Þðt sÞlM ðsÞds dt ð7Þ 0
¼ ðkin lN ÞðtÞ ðkin lN lM ÞðtÞ; with Mð0Þ ¼
Z1
lM ðyÞ
0
Nð0Þ ¼
Z1
Z0 lðxÞ
kin ðsÞdsdx;
ð3Þ
see Appendix 1 for derivation. In case of a constant past of 0 the production kin(t) = kin for t B 0 the initial value (3) results in 0 Nð0Þ ¼ kin T;
ð4Þ
see Appendix 2 for a proof and compare [1]. Summarizing, a distributed lifespan model (DLSM) for a population N(t) with the lifespan PDF lN(s) reads d NðtÞ ¼ kin ðtÞ ðkin lN ÞðtÞ; dt
ð5Þ
with Nð0Þ ¼
0
lN ðxÞ
Z0
kin ðsÞdsdx:
ð6Þ
x
We call (5), (6) the rate of change formulation. Please note that in this formulation for every time point t the convolution over the lifespan interval I has to be calculated. In Fig. 1 (5), (6) is schematically visualized. Additionally, we consider DLSMs with a precursor population (DPLSM), as introduced in [1]. Then (5), (6) serves as a precursor population and individuals expiring
123
ðkin lN ÞðsÞdsdy;
ð8Þ
y
see Appendix 3 for derivation of the initial value. Again for a constant past of the production rate
x
0
Z1
Z0
0 Mð0Þ ¼ kin TM ;
ð9Þ
holds. Please note that M(t) is independent from N(t) and that in (7) the loss term is a double integral. In Fig. 1 (7), (8) is schematically visualized. Standard PKPD software like ADAPT [5] or NONMEM [6] and even general engineering software like MATLAB [7] provide no internal function for the convolution operator. Therefore, the convolution operator has to be approximated by a method provided by the user or an equivalent model formulation is necessary. Krzyzanski et al. [1] presented an approximation of the convolution operator by introducing a large number of additional differential equations. In case of a special variant of the gamma distribution, Nichols et al. [8] constructed an equivalent differential equation system without the convolution but with constraints between the initial values and the model parameters. Freise et al. [2] reported the solution representation of a DLSM without a precursor and applied the less common software WINFUNFIT with additional subroutines. In this paper we explicitly calculate the equivalent integral solution representation of DLSMs without precursor (5), (6) and with precursor (7), (8). The obtained solution could then be simply approximated e.g. by the
J Pharmacokinet Pharmacodyn (2013) 40:639–650
641
Theoretical
Weibull 1 0.9
PDF l(τ) and CDF L(τ)
trapezoidal rule [9]. Therefore, a direct implementation of DLSMs in standard PKPD software is possible. We simulate artificial data to test the approximation of the solution and compare the computational effort with the original rate of change formulation and the approach presented in [1]. Finally, we apply our method to a published model [8] and re-fit hemoglobin (Hgb) data.
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
The main idea in the present paper is to exploit the special structure of Eqs. (5)–(8) and to calculate the solution in an explicit manner. Since (5), (6) do not depend on N(t), the integral solution representation of the rate of change formulation (5), (6) reads NðtÞ ¼
Z1
ð1 LN ðsÞÞkin ðt sÞds;
ð10Þ
0
1 LN ðna hN Þ\a;
ð11Þ
for a given probability significance level a. Condition (11) bounds the probability of large lifespans, see Fig. 2. A reasonable value is a = 0.01. The approximation of (10) reads NðtÞ
Zsend
ð1 LN ðsÞÞkin ðt sÞds ð12Þ
0
hN
na X
ð1 LN ðihN ÞÞkin ðt ihN Þ:
i¼1
For (7), (8) we obtain the integral solution representation MðtÞ ¼
Z1 0
ð1 LM ðxÞÞ
Z1
kin ðt x sÞlN ðsÞdsdx;
ð13Þ
0
see Appendix 6 and e.g. the Riemann sum of (13) reads
4
nh=T
6
Lifespan τ
nα h
8
10
Fig. 2 The solid line is the PDF l(s) and the dashed line is the CDF L(s) of the Weibull distribution with expectation T = 5 and variance r2 = 1. The number of partitions for the lifespan interval [0, T] for a step size h is denoted by n. The minimal number of partitions for the whole lifespan interval up to a significance level a based on condition (11) is denoted by na
0
see Appendix 4 for the proof or the alternative derivation in [2]. The remaining integral over the lifespan interval I in formulation (10) could be simply approximated by standard integration techniques such as the trapezoidal rule (see Appendix 5) or Simpson’s method, see [9]. For illustration purpose, we present the simplest possible approximation (Riemann sum, see [10]) for the realistic lifespan interval eI ¼ ½0; send where send is a natural bound on the lifespan. Let n be a given number of partitions for the interval [0,T] and hN ¼ Tn the step size. The final number of partitions na [ n is the minimal value which fulfills the condition
2
MðtÞ hN hM
nb X na X
ð1 LM ð jhM ÞÞ
j¼1 i¼1
ð14Þ
kin ðt jhM ihN ÞlN ðihN Þ:
Methods In this work we consider three typical examples of continuous distributions, namely Weibull (WBD), gamma (GAD) and lognormal (LND). The WBD is generally applied to describe the lifespan of individuals in statistics [11]. In [12] a combination of two WBDs was presented to describe human mortality and applied by Korell et al. [13] to model red blood cell (RBC) survival. The GAD is e.g. used to describe the lifespan of reticulocytes [1] and also RBCs [8]. The LND is commonly used in clinical PKPD population analysis to describe the distribution of individual model parameter. Each of these distributions consist of two parameters, denoted by a and b, see Table 1. For WBD the gamma function CðsÞ is necessary to calculate the expectation, for LND the Gauss error function erf(s) is needed for the CDF and for GAD the incomplete gamma function c(s) as well as CðsÞ are necessary, see Table 1. These additional functions could e.g. be implemented with the algorithms presented in [14]. ADAPT 5 is a purely Fortran based PKPD software and could simply be extended with additional subroutines and functions provided by the user, see Appendix 7. In [14] Fortran subroutines for cðsÞ; CðsÞ and erf(s) could be found. In ADAPT the tolerance for the simplex algorithm was set to Reqmin = 1E-8 and for the LSODA solver to
123
642
J Pharmacokinet Pharmacodyn (2013) 40:639–650
Table 1 The PDF, CDF, expectation and variance of the Weibull, lognormal and gamma distribution
PDF l(x, a, b) CDF L(x, a, b) Expectation T Variance r2
Weibull
Lognormal
a a x a1 exp bx b b
ðlogðxÞbÞ2 1 pffiffiffiffiffiffiffi exp 2 2 2a x 2pa logðxÞb p ffiffi 0:5 þ 0:5erf 2a 2 exp b þ a2
a 1 exp bx bC 1 þ 1a b2 C 1 þ 2a T 2
Gamma ba a1 CðaÞ x
expðbxÞ
cða; bxÞ CðaÞ a b
a ðexpða2 Þ 1Þ expð2b þ a2 Þ b2 R1 Rx The gamma function is CðxÞ ¼ 0 sx1 expðsÞds; the incomplete gamma function is cðs; xÞ ¼ 0 ys1 expðyÞdy and the Gauss error function Rx 2 2 is erfðxÞ ¼ pffiffip 0 expðs Þds
RTOL = ATOL = 1E-8. In NONMEM 7.2 the $PRED environment was applied and Verbatim Code was used to implement additional functionality. The FOCE method was used for estimation. In both software the trapezoidal rule for the DLSM was implemented. In MATLAB (R2012a) the integral solution representation (10) and the rate of change formulation (5), (6) was implemented. We applied the trapezoidal rule and the MATLAB internal QUAD integrator (adaptive Simpson’s method) to integrate (10). Finally, we solved the rate of change formulation (5), (6) where the convolution was calculated with QUAD and the rate of change was solved with ODE45. The gamma, incomplete gamma and Gauss error function were calculated with the internal commands gamma, gammainc and erf. The optimization function LSQCURVEFIT (Levenberg–Marquardt algorithm) and additionally the optimization method FMINSEARCH (Nelder–Mead simplex) were used, see Example 1 in Results section for more details. The multiple dosing was coded as an explicit function because the drug concentration has to be calculated at different time points in the past, compare (12). In all optimizations an unweighted least squares approach was applied. All source codes are available as supplemental material and at the homepage of the authors. Exemplarily, in Appendix 7 we present the modified ADAPT subroutines and the additional functions of the implementation of the trapezoidal rule for the integral solution representation (10).
Results In this section we demonstrate the power of the presented method. In the first example we produce artificial data in MATLAB and fit those in ADAPT and NONMEM. Additionally, we investigate the computational effort of the presented approximation of the integral solution, the rate of change formulation and the approach presented in [1]. In
123
the second example we apply our method to a model for Hgb measurements in patients with renal disease from [8]. Example 1
PKPD test problem.
To test the presented formulations we construct an artificial PKPD experiment. The PK is described by the Bateman function cðtÞ ¼
doseka ðexpðke tÞ expðka tÞÞ; Vðka ke Þ
ð15Þ
where the parameter are set to ka = 0.7, ke = 0.5 and V = 1. A multiple dosing situation with the time points t = 0 and 7 is considered and four different dosing groups with dose 2 f10; 50; 100; 500g are constructed. The production kin(t) of the population is set in the context of indirect response and reads ( 0 kin for t\0; kin ðtÞ ¼ ð16Þ Smax cðtÞ 0 kin 1 þ SC50 þcðtÞ for t 0: We apply (16) to (4), (5) and obtain the test model d NðtÞ ¼ kin ðtÞ ðkin lN ÞðtÞ; dt
0 Nð0Þ ¼ kin TN ;
ð17Þ
with the model parameter vector 0 h ¼ Smax ; SC50 ; kin ; a; b ; where (a, b) are the distribution related parameter. The WBD is used and the significance level is set to a = 0.01. We investigate three different scenarios. In scenario I, unperturbed data at times 0, 1, 2,…,25 with the parameters h = (0.5, 10, 1, 5.7974, 5.3999) is simulated. In II, this data is moderately perturbed by an additive distributed error with N ð0; 0:04Þ: Finally, we strongly perturbed scenario I by an additive distributed error with N ð0; 0:25Þ; denoted as III. The data is produced with the solution representation (10) where the integrator QUAD is used. The initial estimates h = (0.7, 5, 0.8, 3.88755, 6.6314) are roughly guessed based on the scenario III data. In Fig. 3 the three different produced data scenarios and the corresponding fits are presented. In Table 2 the fitting
J Pharmacokinet Pharmacodyn (2013) 40:639–650
PK
4
Response (Unperturbed) 8 7.5
2
10
7 6.5
N(t)
Drug concentration c(t)
10
0
10
6 5.5
−2
10
5 4.5
−4
10
0
5
10
15
20
4
25
0
5
7.5
7.5
7
7
6.5
6.5
N(t)
8
6
5.5
5
5
4.5
4.5 0
5
10
15
20
15
20
25
20
25
6
5.5
4
10
Response (Strong)
Response (Moderate) 8
N(t)
Fig. 3 Fit of the artificially produced data with the approximation (trapezoidal rule with n = 500) of the solution (10). The upper left panel shows the PK for dose 2 f10; 50; 100; 500g: The upper right panel is the response N(t) of the DLSM fitted to the unperturbed scenario I. The dose = 10 data is denoted by x, the dose = 50 data by squares, the dose = 100 by plus and circles denotes the dose = 500 data. The lower left panel is the fit of the moderately perturbed data II and the lower right panel shows the fit of the strongly perturbed data III
643
4
25
0
5
10
15
Time t
Table 2 Parameter estimates of the approximation (trapezoidal rule with n = 500) of the DLSM for scenarios I (unperturbed), II (moderate perturbed) and III (strong perturbed) obtain from the ADAPT implementation
Parameters
Initial estimates
Scenario I
Scenario II
Scenario III
Model parameters Smax (CV%)
0.7
0.50 (0.32E-2)
0.50 (3.16)
0.494 (7.80)
SC50 (CV%)
5
9.99 (0.11E-1)
9.12 (11.4)
8.01 (28.3)
0 (CV%) kin
0.8
1.00 (0.43E-2)
1.02 (4.21)
1.06 (10.3)
a (CV%)
3.88
5.77 (0.23E-2)
6.36 (3.09)
6.21 (7.01)
b (CV%)
6.63
5.40 (0.63E-2)
5.24 (7.51)
5.04 (17.3)
6 3
5.00 (0.64E-2) 1.01 (0.99E-2)
4.87 (7.65) 0.80 (10.9)
4.68 (17.6) 0.77 (25.6)
Secondary parameters TN (CV%) r2N
(CV%)
Fitting information SofS SofS objective function value, Iter number of iterations of the simplex algorithm
0.2633E-5
3.26162
20.4013
Iter
165
151
100
Time
39 s
33 s
26 s
results of the three different scenarios from the ADAPT implementation of the trapezoidal rule with n = 500 are listed. The NONMEM results are comparable. To get an understanding of the performance of various formulations, three different groups of methods are compared in MATLAB: (i)
Approximation of (10) with the trapezoidal rule or QUAD
(ii)
Approximation of the rate of change formulation (5), (6) (iii) Approximation due to Krzyzanski et al. [1]
The unperturbed scenario I data was used and therefore, in a successful fit the sum of squares is zero. The results are presented in Table 3 in detail. The outcome is that methods of type (i) are extremely superior to those of types (ii) and (iii) in regard of precision and computational effort.
123
644
J Pharmacokinet Pharmacodyn (2013) 40:639–650
Table 3 Numerical performance of the different MATLAB implementations for the scenario I data Formulations
TOL (ODE)
n
Time
Solution Trapeza
100
18 s
7
0.000011
Solution Trapezc
100
5 min 25 s
500d
0.000012
Solution Trapeza
500
1 min 25 s
7
0.000003
25 min
d
0.000003
7
0.000000
Solution Trapez
TOL (QUAD)
c
500
Iter
500
SofS
Solution QUADa
1E-8
37 s
Solution QUADc
1E-8
11 min
500d
1E-8
1E-8
11 min
6
Rate of changea
1E-10
1E-10
3h
27
0.039566
Rate of changea
1E-11
1E-11
16 h
59
0.000000
Rate of change
Rate of change
a
b
0.000000 23.34637
1E-8
1E-8
1 h 11 min
26
0.000153
Rate of changec Method from [1]a
1E-8 1E-8
1E-8
3 h 10 min 1 h 15 min
500d 10
0.000000 2.202447
Method from [1]a
1E-8
100
3 h 10 min
8
2.202672
Method from [1]a
1E-10
50
5 h 40 min
27
1.348316
Method from [1]b
1E-8
50
5 h 30 min
31
0.006383
Method from [1]b
1E-8
100
26 h 30 min
51
0.000000
Method from [1]c
1E-8
50
10 h 30 min
414
0.000000
50
SofS objective function value, Iter number of iterations of the optimization algorithm, TOL (ODE) tolerance used in ODE45 for RTOL = ATOL, TOL (QUAD) tolerance used for the integrator QUAD, n number of partitions of the lifespan interval up to the expectation T in case of the trapezoidal rule or n denotes the number of additional differential equations for the approach from Krzyzanski et al. [1] a
Levenberg–Marquardt with forward difference quotient
b
Levenberg–Marquardt with central difference quotient
c
Simplex algorithm
d
Maximum number of iterations achieved
Additionally, we remark that the use of the GAD or LND leads to similar performance results.
by the GAD (see Table 1) of second order (a = 2). Hence, we have the PDF lM ðsÞ ¼ b2 s expðbsÞ;
Example 2 Stimulation of Hgb by recombinant human erythropoietin (EPO). EPO is a glycoprotein hormone that stimulates RBCs and is primarily produced by kidneys. Hgb is the ironcontaining oxygen-transport metal protein in these RBCs. In renal disease, endogenously produced EPO is insufficient to maintain normal Hgb concentration. Therefore, recombinant human EPO (rHuEPO) is administered to patients to produce the desired concentration of Hgb. Nichols et al. [8] presented a rich data set of Hgb measurements from renal disease patients over 460 days with about three measurements a week. At every measurement time point rHuEPO was administered to patients. The amount of administered rHuEPO is known but drug concentration profiles were not available. For this example we digitized the data from patient 1. The authors modeled the data by a DLSM with precursor of the form (5)–(8) where the precursor N(t) describes hematopoietic stem cells which maturate after a constant lifespan TN into RBCs denoted by M(t). The lifespan of the RBCs was described
123
ð18Þ
and the CDF could be explicitly calculated and reads LM ðsÞ ¼ 1 ð1 þ bsÞ expðsbÞ:
ð19Þ
EPO concentration in plasma is the sum of endogenous EPO Een, a constant value, and exogenous rHuEPO E(t). The PK of rHuEPO is described by nonlinear elimination d Vmax EðtÞ EðtÞ ¼ ; Eð0Þ ¼ 0; dt Km þ EðtÞ and the total concentration of EPO reads cðtÞ ¼ EðtÞ þ Een : The stimulatory effect is kin ðtÞ ¼
Smax cðtÞ SC50 þ cðtÞ
for
t 0;
Een 0 and kin ðtÞ ¼ kin ¼ SCSmax for t \ 0. The final model from 50 þEen [8] reads
d MðtÞ ¼ kin ðt TN Þ ðkin lM Þðt TN Þ; dt
ð20Þ
J Pharmacokinet Pharmacodyn (2013) 40:639–650
645
Hemoglobin H(t) (g/dL)
15
Table 4 Parameter estimates of the approximation (trapezoidal rule [n = 500]) of the hemoglobin model (20)
14
Parameters
Units
Estimates
Vmax
IU/day
6,890a
Km
IU
1,042a
13
Model parameters 12
11
10
9 0
50
100
150
200
250
300
350
400
450
500
Time t (day)
Fig. 4 Fit of the hemoglobin data with the approximation (trapezoidal rule with n = 500) of the solution (22) of the distributed lifespan model with precursor (20)
HðtÞ ¼ KH MðtÞ;
ð21Þ
where M(t) describes the RBC and H(t) the Hgb concentration. The value KH = 29.5 g/dL is the average amount of Hgb per RBC. The model parameters are h ¼ ðVmax ; Km ; Hen ; Smax ; SC50 ; TN ; bÞ: where Hen is the endogenous level of Hgb due to Een, see [8]. In [8] a three dimensional ordinary differential equation system for (20) based on (18) was constructed with constraints between initial values and model parameter. This reformulation was applied in [8] to fit the Hgb data. Unfortunately, in their reformulation of (20) the initial value is only formulated as a constraint between the initial values of the additional ODEs and not explicitly given. In DLSMs with known inflow for t B 0, according to (7)–(9) there is no freedom to choose 0 M(0). Therefore, due to (9) we have M(0) = kin TM. In this work we apply the solution representation (13) of the precursor model (7), (8) to (20) and obtain Z1 MðtÞ ¼ ð1 LM ðsÞÞkin ðt TN sÞds; ð22Þ 0
where LM(s) is given by (19). The solution (22) is approximated by the trapezoidal rule with n = 500 and a = 0.01. In addition, we apply the PK averaging technique presented in [8] because c(t) has an impulse-like behavior (186 doses over 460 days with rapid clearance). Therefore, Nichols et al. suggested to average the concentration in certain time intervals by Smax h Km Een kin ðdi Þ ¼ Km SC50 Ti Vmax SC50 þ Een SC50 þ Een þ di ln SC50 þ Een i Een Smax Een þ 1 ; ð23Þ di þ SC50 þ Een SC50 þ Een
Hen (CV%)
g/dL
10.26 (1.97)
Smax
cell/day
0.207E-1a
SC50 (CV%)
IU
86.71 (83.1)
TN (CV%) b (CV%)
day
6.98 (10.0) 0.605E-1 (10.0)
TM (CV%)
day
33.1 (10.0)
r2M
day2
547 (20.0)
Secondary parameters (CV%)
Fitting information
a
SofS
68.84
Iter
129
fixed to the values presented in [8]
where Ti is a time interval and di the dose at the time point ti. However, in the following some discussion is necessary. In [8] it is not obvious how Ti was set in (23). Here we choose the interval length Ti for a given dose di as the difference between the time point of the next dose di?1 and the time point of di. Also the PK parameters (Vmax,Km) (no PK data available) are estimated just from the Hgb PD data. This could lead to an over-parameterization of the model. Therefore, we fix Vmax, Km and also Smax to the values presented in the original work. In Fig. 4 the fit of the approximation of (22) with the trapezoidal rule (n = 500) is presented and in Table 4 the parameter estimates of (SC50, Hen, TN, b) are shown from the ADAPT implementation. The fit and the estimates do not exactly coincide due to the differences with the setting of the initial value M(0). Unfortunately, in [8] no derivation of their reformulation of (20) is presented and no information about the optimization process is given.
Discussion Previously published population models with a distributed lifespan (see [1]) based on the convolution operator were mathematically observed. The appearance of the convolution in the rate of change formulation makes an implementation extremely difficult in standard PKPD software. In the presented work we calculated the integral solution representation of DLSMs with and without a precursor population for arbitrary distributions and approximated these representations e.g. by the trapezoidal rule. Therefore, these models could now be simply implemented
123
646
J Pharmacokinet Pharmacodyn (2013) 40:639–650
in every PKPD software. In addition, the presented results are formulated in the general case of a non-constant past for the production kin(t). Further we investigated computational properties of different model formulations. We demonstrated that approximating the calculated solution representation instead of solving the rate of change formulation or the approximation presented in [1], results in a tremendous computational time advantage. Moreover, it was observed that our approach avoids numerical optimization problems. Finally, we applied our method to a published Hgb model and data from [8]. For the sake of completeness we mention that if the production kin(t) is not explicitly given, the advantage of our method will decrease because an additional ODE has to be simultaneously solved. This leads to a relevant increase in computational time. To overcome this issue other approximation techniques are necessary. For the special case of lifespan models with a point distribution it is already proven that an approximation with transit compartments exists, see [15, 16]. However, for lifespan models with an arbitrary continuous distribution an approximation by transit compartments is topic of our next research.
Then lim NðtÞ ¼ 0
t!1
, Nð0Þ ¼
Z1
lN ðxÞ
Z1
0
kf in ðs xÞdsdx:
0
Using (24) Nð0Þ ¼
¼
Z1
lN ðxÞ
Zx
0
0
Z1
Z0
lN ðxÞ
kin ðs xÞdsdx
kin ðsÞdsdx:
x
0
Appendix 2 0 We consider a constant past kin(t) = kin for t B 0 and obtain for the initial values
Nð0Þ ¼
Acknowledgments The present project is supported by the National Research Fund, Luxembourg, and cofunded under the Marie Curie Actions of the European Commission (FP7-COFUND).
Z1
Z0
lN ðxÞ
0 kin dsdx
x
0 0 ¼ kin
Z1
0 lN ðxÞxdx ¼ kin TN ;
0
Appendices
and
Appendix 1 Mð0Þ ¼ To estimate the amount of individuals at the initial population state N(0) for the DLSM (5), (6) we define the production term kin ðtÞ t 0; f kin ðtÞ ¼ ð24Þ 0 t [ 0: Then the appropriate initial value N(0) is determined by the claim limt!1 NðtÞ ¼ 0: Substituting (24) into (5) results in d NðtÞ ¼ kf in ðtÞ dt
Z1
kf in ðt xÞlN ðxÞdx:
ð25Þ
¼
¼
Z1
lM ðyÞ
Z0
0
y
Z1
Z0
lM ðyÞ
0
y
Z1
Z0
lM ðyÞ
0 kin lN ðsÞdsdy
0 kin
Z1
lN ðxÞdxdsdy
0 0 0 kin dsdy ¼ kin TM :
y
0
Appendix 3
0
To estimate the amount M(0) of the follow-up population for the DPLSM (7), (8) a similar approach as in Appendix 1 is applied. We define the production
The general solution of (25) reads NðtÞ ¼ Nð0Þ þ
¼ Nð0Þ
Zt 0 Z1 0
123
kf in ðsÞds
Z t Z1 0
lN ðxÞ
Zt 0
kf in ðs xÞlN ðxÞdxds
0
kf in ðs xÞdsdx:
kf in ðtÞ ¼
ðkin lÞðtÞ 0
t 0; t [ 0;
ð26Þ
and use the condition limt!1 MðtÞ ¼ 0 to determine M(0). Substituting (26) into (7) results in
J Pharmacokinet Pharmacodyn (2013) 40:639–650
MðtÞ ¼ Mð0Þ þ
¼ Mð0Þ
Zt
kf in ðsÞds
0 Z1
Z t Z1 0
lM ðyÞ
0
Zt
647
Appendix 5 kf in ðs yÞlM ðyÞdyds Generally, for an arbitrary continuous function f, the trapezoidal rule reads
0
kf in ðs yÞdsdy;
Zsend
0
f ðsÞds
s0
and therefore, Mð0Þ ¼
¼
Z1
lM ðyÞ
Z1
0
0
Z1
Z0
lM ðyÞ
with the equidistant partition sj = s0 ? (j - 1)h with h ¼ for j = 1,…,m ? 1 of the lifespan interval [s0, send].
kf in ðs yÞdsdy ð27Þ ðkin lN ÞðsÞdsdy:
Appendix 4
Appendix 6 Differentiation of (13) with respect to time t gives Z1 Z1 d 0 MðtÞ ¼ ð1 LM ðxÞÞ kin ðt x sÞlN ðsÞdsdx dt 0
Differentiation of (10) with respect to time t gives d NðtÞ ¼ dt
Z1
¼
0
Z1
lN ðsÞ
0
ð1
0 LN ðsÞÞkin ðt
ð28Þ
send s0 m
y
0
m hX ð f ðsiþ1 Þ þ f ðsi ÞÞ; 2 i¼1
Z1
0
ð1 LM ðxÞÞkin ðt x sÞdxds:
0
ð29Þ
sÞds
0
¼ ½ð1 LN ðsÞÞkin ðt sÞs¼1 s¼0 Z1 ðlN ðsÞÞðkin ðt sÞÞds
The inner integral in (29) could be written as Z1
0
ð1 LM ðxÞÞkin ðt x sÞdx
0 0
¼ kin ðtÞ
Z1
lN ðsÞkin ðt sÞds
0
With (10) we obtain for the initial value Nð0Þ ¼
2 ¼ 4ð1 LN ðsÞÞ
Zs
Z1
ðlN ðsÞÞ
0
¼
0
kin ðt x sÞlM ðxÞdx:
ð31Þ
lN ðsÞ
Zs
Substitution of (31) into (29) gives 3s¼1
kin ðxÞdx5
0
Z1
¼ kin ðt sÞ
Z1 0
ð1 LN ðsÞÞkin ðsÞds
0
ð30Þ
0
¼ kin ðtÞ ðkin lN ÞðtÞ:
Z1
¼ ½ð1 LM ðxÞÞkin ðt x sÞx¼1 x¼0 Z1 kin ðt x sÞlM ðxÞdx
d MðtÞ ¼ dt
Z1 0
s¼0
kin ðxÞdxds
lN ðsÞkin ðt sÞds
Z1 0
lN ðsÞ
Z1
kin ðt x sÞlM ðxÞdxds
0
0
Z0 s
¼ ðkin lN ÞðtÞ kin ðxÞdxds:
Z1
lN ðsÞððkin lM Þðt sÞÞds
0
¼ ðkin lN ÞðtÞ ðkin lN lM ÞðtÞ:
123
648
J Pharmacokinet Pharmacodyn (2013) 40:639–650
For the initial value we obtain Mð0Þ ¼
Z1
ð1 LM ðxÞÞ
Z1
0
Z x Z1
¼ 4ð1 LM ðxÞÞ
0
þ
Z1
lM ðxÞ
0
¼
¼
kin ðx sÞlN ðsÞdsdx
0
2
Z1
Z x Z1 0
lM ðxÞ
Zx
kin ðy sÞlN ðsÞdsdy5
0
x¼0
kin ðy sÞlN ðsÞdsdydx
ðkin lN ÞðyÞdydx
0
0
Z0
0
3x¼1
0
Z1
lM ðxÞ
CDF L(x) of the WBD:
ðkin lN ÞðyÞdydx:
x
Appendix 7 In the following we present the implementation of the trapezoidal rule for (10). Only additional functions and modified ADAPT subroutines are presented. For the complete code see the supplemental material. Logarithm of the gamma function CðxÞ based on [14]:
123
Multiple dosing function for c(t):
J Pharmacokinet Pharmacodyn (2013) 40:639–650
Production kin(t):
649
The ADAPT subroutine OUTPUT calls the trapezoidal rule for the different dosing groups.
Integrand in (10):
Trapezoidal rule:
123
650
References 1. Krzyzanski W, Woo S, Jusko WJ (2006) Pharmacodynamic models for agents that alter production of natural cells with various distributions of lifespans. J Pharmacokinet Pharmacodyn 33(2):125–166
123
J Pharmacokinet Pharmacodyn (2013) 40:639–650 2. Freise K, Widness J, Schmidt R, Veng-Pedersen P (2008) Modeling time variant distributions of cellular lifespans: increases in circulating reticulocyte lifespans following double phlebotomies in sheep. J Pharmacokinet Pharmacodyn 35(3):285–323 3. Ramakrishnan R, Cheung WK, Wacholtz MC, Minton N, Jusko WJ (2004) Pharmacokinetic and pharmacodynamic modeling of recombinant human erythropoietin after single and multiple doses in healthy volunteers. J Clin Pharmacol 44:991–1002 4. Ait-Oudhia S, Scherrmann JM, Krzyzanski W (2010) Simultaneous pharmacokinetics/pharmacodynamics modeling of recombinant human erythropoietin upon multiple intravenous dosing in rats. J Pharmacol Exp Ther 334:897–910 5. D’Argenio DZ, Schumitzky A, Wang X (2009) ADAPT 5 user’s guide: pharmacokinetic/pharmacodynamic systems analysis software. Biomedical Simulations Resource, Los Angeles 6. Beal S, Sheiner LB, Boeckmann A, Bauer RJ (2009) NONMEM user’s guides. Icon Development Solutions, Ellicott City 7. MATLAB Release (2012) The MathWorks, Inc., Natick 8. Nichols B, Shrestha RP, Horowitz J, Hollot CV, Germain MJ, Gaweda AE, Chait Y (2011) Simplification of an erythropoiesis model for design of anemia management protocols in end stage renal disease. In: 33rd Annual international conference of the IEEE EMBS Boston, Massachusetts, USA 9. Atkinson K (1989) An introduction to numerical analysis, 2nd edn. Wiley, New York 10. Thomas GB, Finney RL (1996) Calculus and analytic geometry, 9th edn. Addison-Wesley, Reading 11. Marshall AW, Olkin I (2007) Life distributions. Springer, New York 12. Bebbinton M, Lai CD, Zitikis R (2007) Modeling human mortality using mixtures of bathtub shaped failure distributions. J Theor Biol 245(3):528–538 13. Korell J, Coulter CV, Duffull SB (2011) A statistical model for red blood cell survival. J Theor Biol 268(1):39–49 14. Press WH, Teukolsky SA, Vetterling WT, Flannery BP (1992) Numerical recipes in Fortran 77, the art of scientific computing, 2nd edn. Cambridge University Press, New York 15. Krzyzanski W (2011) Interpretation of transit compartments pharmacodynamic models as lifespan based indirect response models. J Pharmacokinet Pharmacodyn 38(2):179–204 16. Koch G, Schropp J (2012) General relationship between transit compartments and lifespan models. J Pharmacokinet Pharmacodyn 39(4):343–355
ORIGINAL PAPER
Solution and implementation of distributed lifespan models Gilbert Koch • Johannes Schropp
Received: 25 July 2013 / Accepted: 8 October 2013 / Published online: 1 November 2013 Springer Science+Business Media New York 2013
Abstract We consider a population where every individual has a unique lifespan. After expiring of its lifespan the individual has to leave the population. A realistic approach to describe these lifespans is by a continuous distribution. Such distributed lifespan models (DLSMs) were introduced earlier in the indirect response context and consist of the mathematical convolution operator to describe the rate of change. Therefore, DLSMs could not directly be implemented in standard PKPD software. In this work we present the solution representation of DLSMs with and without a precursor population and an implementation strategy for DLSMs in ADAPT, NONMEM and MATLAB. We fit hemoglobin measurements from literature and investigate computational properties. Keywords Lifespan Distribution Convolution Population
After expiration of its lifespan the individual has to leave the population and maybe maturate into another population. A realistic approach to describe lifespans is by a continuous distribution. Such models were introduced by Krzyzanski et al. [1] to PKPD modeling in the indirect response context in 2006. In [1, 2] these models where applied to describe the stimulation of reticulocytes. In [3, 4] simpler models where every individual has the same lifespan are used to describe erythropoiesis. Let T be a random variable describing the distribution of the lifespan s of individuals based on a probability density function (PDF) l(s). The only assumption is that negative lifespans are not possible, more precisely l(s) = 0 for s \ 0. The expectation of the lifespan distribution is denoted by Z1 T ¼ E½T ¼ slðsÞds; 0
Introduction
the variance is
We consider the development of a population N(t) in time t where every individual has its own and unique lifespan s.
r2 ¼ Var½T ¼
Electronic supplementary material The online version of this article (doi:10.1007/s10928-013-9336-y) contains supplementary material, which is available to authorized users. G. Koch (&) Department of Pharmaceutical Sciences, State University of New York at Buffalo, 403 Kapoor Hall, Buffalo, NY 14214, USA e-mail: [email protected] J. Schropp Department of Mathematics and Statistics, Universita¨t Konstanz, PO Box 195, 78457 Konstanz, Germany
Z1
s2 lðsÞds T 2 ;
0
and the cumulative distribution function (CDF) is defined by LðsÞ ¼
Zs lðsÞds for
s 0:
0
In the special case of a point distribution (i.e. every individual has the same lifespan T) we have l(s) = d(s - T) where d is the Dirac delta function. A model describing the rate of change of a population reads
123
640
J Pharmacokinet Pharmacodyn (2013) 40:639–650
d NðtÞ ¼ kin ðtÞ kout ðtÞ; dt
Nð0Þ ¼ N 0 ;
ð1Þ
where the production (birth) rate is described by kin(t) and kout(t) denotes the loss (death) rate. The initial amount of individuals in the population at time point t = 0 is denoted by N0. Please note that only zero order terms appear in (1). The lifespan distribution is included into (1) via the convolution operator Z1 kout ðtÞ ¼ ðkin lÞðtÞ ¼ kin ðt sÞlðsÞds; ð2Þ
Fig. 1 The upper scheme is the basic distributed lifespan model, DLSM (5), (6) and the lower scheme shows the distributed lifespan model with precursor, DPLSM (7), (8)
0
see [1] for derivation. Hence, to calculate the outflow kout(t) at a time point t an integral over the lifespan interval I ¼ ½0; 1 has to be solved. Note that in (2) the past t \ 0 of the production kin(t) plays an important role and has to be defined. In this work we admit a non-constant past for the production kin(t) for t \ 0 which is an extension of traditional PKPD modeling and covers more pharmacological situations. Please note that (2) reduces to kout(t) = kin(t - T) in case of a point distributed lifespan, see [1]. The initial value N(0) in (1) describes the amount of born and already died individuals in the past up to the starting time point t = 0 and is calculated by
from that population maturate into a follow-up population M(t). Mathematically, this means that kin and lN from (5), (6) have to be replaced by kin * lN and lM. Hence, for a DPLSM we additionally obtain Z1 d MðtÞ ¼ ðkin lN ÞðtÞ ðkin lN Þðt sÞlM ðsÞds dt ð7Þ 0
¼ ðkin lN ÞðtÞ ðkin lN lM ÞðtÞ; with Mð0Þ ¼
Z1
lM ðyÞ
0
Nð0Þ ¼
Z1
Z0 lðxÞ
kin ðsÞdsdx;
ð3Þ
see Appendix 1 for derivation. In case of a constant past of 0 the production kin(t) = kin for t B 0 the initial value (3) results in 0 Nð0Þ ¼ kin T;
ð4Þ
see Appendix 2 for a proof and compare [1]. Summarizing, a distributed lifespan model (DLSM) for a population N(t) with the lifespan PDF lN(s) reads d NðtÞ ¼ kin ðtÞ ðkin lN ÞðtÞ; dt
ð5Þ
with Nð0Þ ¼
0
lN ðxÞ
Z0
kin ðsÞdsdx:
ð6Þ
x
We call (5), (6) the rate of change formulation. Please note that in this formulation for every time point t the convolution over the lifespan interval I has to be calculated. In Fig. 1 (5), (6) is schematically visualized. Additionally, we consider DLSMs with a precursor population (DPLSM), as introduced in [1]. Then (5), (6) serves as a precursor population and individuals expiring
123
ðkin lN ÞðsÞdsdy;
ð8Þ
y
see Appendix 3 for derivation of the initial value. Again for a constant past of the production rate
x
0
Z1
Z0
0 Mð0Þ ¼ kin TM ;
ð9Þ
holds. Please note that M(t) is independent from N(t) and that in (7) the loss term is a double integral. In Fig. 1 (7), (8) is schematically visualized. Standard PKPD software like ADAPT [5] or NONMEM [6] and even general engineering software like MATLAB [7] provide no internal function for the convolution operator. Therefore, the convolution operator has to be approximated by a method provided by the user or an equivalent model formulation is necessary. Krzyzanski et al. [1] presented an approximation of the convolution operator by introducing a large number of additional differential equations. In case of a special variant of the gamma distribution, Nichols et al. [8] constructed an equivalent differential equation system without the convolution but with constraints between the initial values and the model parameters. Freise et al. [2] reported the solution representation of a DLSM without a precursor and applied the less common software WINFUNFIT with additional subroutines. In this paper we explicitly calculate the equivalent integral solution representation of DLSMs without precursor (5), (6) and with precursor (7), (8). The obtained solution could then be simply approximated e.g. by the
J Pharmacokinet Pharmacodyn (2013) 40:639–650
641
Theoretical
Weibull 1 0.9
PDF l(τ) and CDF L(τ)
trapezoidal rule [9]. Therefore, a direct implementation of DLSMs in standard PKPD software is possible. We simulate artificial data to test the approximation of the solution and compare the computational effort with the original rate of change formulation and the approach presented in [1]. Finally, we apply our method to a published model [8] and re-fit hemoglobin (Hgb) data.
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
The main idea in the present paper is to exploit the special structure of Eqs. (5)–(8) and to calculate the solution in an explicit manner. Since (5), (6) do not depend on N(t), the integral solution representation of the rate of change formulation (5), (6) reads NðtÞ ¼
Z1
ð1 LN ðsÞÞkin ðt sÞds;
ð10Þ
0
1 LN ðna hN Þ\a;
ð11Þ
for a given probability significance level a. Condition (11) bounds the probability of large lifespans, see Fig. 2. A reasonable value is a = 0.01. The approximation of (10) reads NðtÞ
Zsend
ð1 LN ðsÞÞkin ðt sÞds ð12Þ
0
hN
na X
ð1 LN ðihN ÞÞkin ðt ihN Þ:
i¼1
For (7), (8) we obtain the integral solution representation MðtÞ ¼
Z1 0
ð1 LM ðxÞÞ
Z1
kin ðt x sÞlN ðsÞdsdx;
ð13Þ
0
see Appendix 6 and e.g. the Riemann sum of (13) reads
4
nh=T
6
Lifespan τ
nα h
8
10
Fig. 2 The solid line is the PDF l(s) and the dashed line is the CDF L(s) of the Weibull distribution with expectation T = 5 and variance r2 = 1. The number of partitions for the lifespan interval [0, T] for a step size h is denoted by n. The minimal number of partitions for the whole lifespan interval up to a significance level a based on condition (11) is denoted by na
0
see Appendix 4 for the proof or the alternative derivation in [2]. The remaining integral over the lifespan interval I in formulation (10) could be simply approximated by standard integration techniques such as the trapezoidal rule (see Appendix 5) or Simpson’s method, see [9]. For illustration purpose, we present the simplest possible approximation (Riemann sum, see [10]) for the realistic lifespan interval eI ¼ ½0; send where send is a natural bound on the lifespan. Let n be a given number of partitions for the interval [0,T] and hN ¼ Tn the step size. The final number of partitions na [ n is the minimal value which fulfills the condition
2
MðtÞ hN hM
nb X na X
ð1 LM ð jhM ÞÞ
j¼1 i¼1
ð14Þ
kin ðt jhM ihN ÞlN ðihN Þ:
Methods In this work we consider three typical examples of continuous distributions, namely Weibull (WBD), gamma (GAD) and lognormal (LND). The WBD is generally applied to describe the lifespan of individuals in statistics [11]. In [12] a combination of two WBDs was presented to describe human mortality and applied by Korell et al. [13] to model red blood cell (RBC) survival. The GAD is e.g. used to describe the lifespan of reticulocytes [1] and also RBCs [8]. The LND is commonly used in clinical PKPD population analysis to describe the distribution of individual model parameter. Each of these distributions consist of two parameters, denoted by a and b, see Table 1. For WBD the gamma function CðsÞ is necessary to calculate the expectation, for LND the Gauss error function erf(s) is needed for the CDF and for GAD the incomplete gamma function c(s) as well as CðsÞ are necessary, see Table 1. These additional functions could e.g. be implemented with the algorithms presented in [14]. ADAPT 5 is a purely Fortran based PKPD software and could simply be extended with additional subroutines and functions provided by the user, see Appendix 7. In [14] Fortran subroutines for cðsÞ; CðsÞ and erf(s) could be found. In ADAPT the tolerance for the simplex algorithm was set to Reqmin = 1E-8 and for the LSODA solver to
123
642
J Pharmacokinet Pharmacodyn (2013) 40:639–650
Table 1 The PDF, CDF, expectation and variance of the Weibull, lognormal and gamma distribution
PDF l(x, a, b) CDF L(x, a, b) Expectation T Variance r2
Weibull
Lognormal
a a x a1 exp bx b b
ðlogðxÞbÞ2 1 pffiffiffiffiffiffiffi exp 2 2 2a x 2pa logðxÞb p ffiffi 0:5 þ 0:5erf 2a 2 exp b þ a2
a 1 exp bx bC 1 þ 1a b2 C 1 þ 2a T 2
Gamma ba a1 CðaÞ x
expðbxÞ
cða; bxÞ CðaÞ a b
a ðexpða2 Þ 1Þ expð2b þ a2 Þ b2 R1 Rx The gamma function is CðxÞ ¼ 0 sx1 expðsÞds; the incomplete gamma function is cðs; xÞ ¼ 0 ys1 expðyÞdy and the Gauss error function Rx 2 2 is erfðxÞ ¼ pffiffip 0 expðs Þds
RTOL = ATOL = 1E-8. In NONMEM 7.2 the $PRED environment was applied and Verbatim Code was used to implement additional functionality. The FOCE method was used for estimation. In both software the trapezoidal rule for the DLSM was implemented. In MATLAB (R2012a) the integral solution representation (10) and the rate of change formulation (5), (6) was implemented. We applied the trapezoidal rule and the MATLAB internal QUAD integrator (adaptive Simpson’s method) to integrate (10). Finally, we solved the rate of change formulation (5), (6) where the convolution was calculated with QUAD and the rate of change was solved with ODE45. The gamma, incomplete gamma and Gauss error function were calculated with the internal commands gamma, gammainc and erf. The optimization function LSQCURVEFIT (Levenberg–Marquardt algorithm) and additionally the optimization method FMINSEARCH (Nelder–Mead simplex) were used, see Example 1 in Results section for more details. The multiple dosing was coded as an explicit function because the drug concentration has to be calculated at different time points in the past, compare (12). In all optimizations an unweighted least squares approach was applied. All source codes are available as supplemental material and at the homepage of the authors. Exemplarily, in Appendix 7 we present the modified ADAPT subroutines and the additional functions of the implementation of the trapezoidal rule for the integral solution representation (10).
Results In this section we demonstrate the power of the presented method. In the first example we produce artificial data in MATLAB and fit those in ADAPT and NONMEM. Additionally, we investigate the computational effort of the presented approximation of the integral solution, the rate of change formulation and the approach presented in [1]. In
123
the second example we apply our method to a model for Hgb measurements in patients with renal disease from [8]. Example 1
PKPD test problem.
To test the presented formulations we construct an artificial PKPD experiment. The PK is described by the Bateman function cðtÞ ¼
doseka ðexpðke tÞ expðka tÞÞ; Vðka ke Þ
ð15Þ
where the parameter are set to ka = 0.7, ke = 0.5 and V = 1. A multiple dosing situation with the time points t = 0 and 7 is considered and four different dosing groups with dose 2 f10; 50; 100; 500g are constructed. The production kin(t) of the population is set in the context of indirect response and reads ( 0 kin for t\0; kin ðtÞ ¼ ð16Þ Smax cðtÞ 0 kin 1 þ SC50 þcðtÞ for t 0: We apply (16) to (4), (5) and obtain the test model d NðtÞ ¼ kin ðtÞ ðkin lN ÞðtÞ; dt
0 Nð0Þ ¼ kin TN ;
ð17Þ
with the model parameter vector 0 h ¼ Smax ; SC50 ; kin ; a; b ; where (a, b) are the distribution related parameter. The WBD is used and the significance level is set to a = 0.01. We investigate three different scenarios. In scenario I, unperturbed data at times 0, 1, 2,…,25 with the parameters h = (0.5, 10, 1, 5.7974, 5.3999) is simulated. In II, this data is moderately perturbed by an additive distributed error with N ð0; 0:04Þ: Finally, we strongly perturbed scenario I by an additive distributed error with N ð0; 0:25Þ; denoted as III. The data is produced with the solution representation (10) where the integrator QUAD is used. The initial estimates h = (0.7, 5, 0.8, 3.88755, 6.6314) are roughly guessed based on the scenario III data. In Fig. 3 the three different produced data scenarios and the corresponding fits are presented. In Table 2 the fitting
J Pharmacokinet Pharmacodyn (2013) 40:639–650
PK
4
Response (Unperturbed) 8 7.5
2
10
7 6.5
N(t)
Drug concentration c(t)
10
0
10
6 5.5
−2
10
5 4.5
−4
10
0
5
10
15
20
4
25
0
5
7.5
7.5
7
7
6.5
6.5
N(t)
8
6
5.5
5
5
4.5
4.5 0
5
10
15
20
15
20
25
20
25
6
5.5
4
10
Response (Strong)
Response (Moderate) 8
N(t)
Fig. 3 Fit of the artificially produced data with the approximation (trapezoidal rule with n = 500) of the solution (10). The upper left panel shows the PK for dose 2 f10; 50; 100; 500g: The upper right panel is the response N(t) of the DLSM fitted to the unperturbed scenario I. The dose = 10 data is denoted by x, the dose = 50 data by squares, the dose = 100 by plus and circles denotes the dose = 500 data. The lower left panel is the fit of the moderately perturbed data II and the lower right panel shows the fit of the strongly perturbed data III
643
4
25
0
5
10
15
Time t
Table 2 Parameter estimates of the approximation (trapezoidal rule with n = 500) of the DLSM for scenarios I (unperturbed), II (moderate perturbed) and III (strong perturbed) obtain from the ADAPT implementation
Parameters
Initial estimates
Scenario I
Scenario II
Scenario III
Model parameters Smax (CV%)
0.7
0.50 (0.32E-2)
0.50 (3.16)
0.494 (7.80)
SC50 (CV%)
5
9.99 (0.11E-1)
9.12 (11.4)
8.01 (28.3)
0 (CV%) kin
0.8
1.00 (0.43E-2)
1.02 (4.21)
1.06 (10.3)
a (CV%)
3.88
5.77 (0.23E-2)
6.36 (3.09)
6.21 (7.01)
b (CV%)
6.63
5.40 (0.63E-2)
5.24 (7.51)
5.04 (17.3)
6 3
5.00 (0.64E-2) 1.01 (0.99E-2)
4.87 (7.65) 0.80 (10.9)
4.68 (17.6) 0.77 (25.6)
Secondary parameters TN (CV%) r2N
(CV%)
Fitting information SofS SofS objective function value, Iter number of iterations of the simplex algorithm
0.2633E-5
3.26162
20.4013
Iter
165
151
100
Time
39 s
33 s
26 s
results of the three different scenarios from the ADAPT implementation of the trapezoidal rule with n = 500 are listed. The NONMEM results are comparable. To get an understanding of the performance of various formulations, three different groups of methods are compared in MATLAB: (i)
Approximation of (10) with the trapezoidal rule or QUAD
(ii)
Approximation of the rate of change formulation (5), (6) (iii) Approximation due to Krzyzanski et al. [1]
The unperturbed scenario I data was used and therefore, in a successful fit the sum of squares is zero. The results are presented in Table 3 in detail. The outcome is that methods of type (i) are extremely superior to those of types (ii) and (iii) in regard of precision and computational effort.
123
644
J Pharmacokinet Pharmacodyn (2013) 40:639–650
Table 3 Numerical performance of the different MATLAB implementations for the scenario I data Formulations
TOL (ODE)
n
Time
Solution Trapeza
100
18 s
7
0.000011
Solution Trapezc
100
5 min 25 s
500d
0.000012
Solution Trapeza
500
1 min 25 s
7
0.000003
25 min
d
0.000003
7
0.000000
Solution Trapez
TOL (QUAD)
c
500
Iter
500
SofS
Solution QUADa
1E-8
37 s
Solution QUADc
1E-8
11 min
500d
1E-8
1E-8
11 min
6
Rate of changea
1E-10
1E-10
3h
27
0.039566
Rate of changea
1E-11
1E-11
16 h
59
0.000000
Rate of change
Rate of change
a
b
0.000000 23.34637
1E-8
1E-8
1 h 11 min
26
0.000153
Rate of changec Method from [1]a
1E-8 1E-8
1E-8
3 h 10 min 1 h 15 min
500d 10
0.000000 2.202447
Method from [1]a
1E-8
100
3 h 10 min
8
2.202672
Method from [1]a
1E-10
50
5 h 40 min
27
1.348316
Method from [1]b
1E-8
50
5 h 30 min
31
0.006383
Method from [1]b
1E-8
100
26 h 30 min
51
0.000000
Method from [1]c
1E-8
50
10 h 30 min
414
0.000000
50
SofS objective function value, Iter number of iterations of the optimization algorithm, TOL (ODE) tolerance used in ODE45 for RTOL = ATOL, TOL (QUAD) tolerance used for the integrator QUAD, n number of partitions of the lifespan interval up to the expectation T in case of the trapezoidal rule or n denotes the number of additional differential equations for the approach from Krzyzanski et al. [1] a
Levenberg–Marquardt with forward difference quotient
b
Levenberg–Marquardt with central difference quotient
c
Simplex algorithm
d
Maximum number of iterations achieved
Additionally, we remark that the use of the GAD or LND leads to similar performance results.
by the GAD (see Table 1) of second order (a = 2). Hence, we have the PDF lM ðsÞ ¼ b2 s expðbsÞ;
Example 2 Stimulation of Hgb by recombinant human erythropoietin (EPO). EPO is a glycoprotein hormone that stimulates RBCs and is primarily produced by kidneys. Hgb is the ironcontaining oxygen-transport metal protein in these RBCs. In renal disease, endogenously produced EPO is insufficient to maintain normal Hgb concentration. Therefore, recombinant human EPO (rHuEPO) is administered to patients to produce the desired concentration of Hgb. Nichols et al. [8] presented a rich data set of Hgb measurements from renal disease patients over 460 days with about three measurements a week. At every measurement time point rHuEPO was administered to patients. The amount of administered rHuEPO is known but drug concentration profiles were not available. For this example we digitized the data from patient 1. The authors modeled the data by a DLSM with precursor of the form (5)–(8) where the precursor N(t) describes hematopoietic stem cells which maturate after a constant lifespan TN into RBCs denoted by M(t). The lifespan of the RBCs was described
123
ð18Þ
and the CDF could be explicitly calculated and reads LM ðsÞ ¼ 1 ð1 þ bsÞ expðsbÞ:
ð19Þ
EPO concentration in plasma is the sum of endogenous EPO Een, a constant value, and exogenous rHuEPO E(t). The PK of rHuEPO is described by nonlinear elimination d Vmax EðtÞ EðtÞ ¼ ; Eð0Þ ¼ 0; dt Km þ EðtÞ and the total concentration of EPO reads cðtÞ ¼ EðtÞ þ Een : The stimulatory effect is kin ðtÞ ¼
Smax cðtÞ SC50 þ cðtÞ
for
t 0;
Een 0 and kin ðtÞ ¼ kin ¼ SCSmax for t \ 0. The final model from 50 þEen [8] reads
d MðtÞ ¼ kin ðt TN Þ ðkin lM Þðt TN Þ; dt
ð20Þ
J Pharmacokinet Pharmacodyn (2013) 40:639–650
645
Hemoglobin H(t) (g/dL)
15
Table 4 Parameter estimates of the approximation (trapezoidal rule [n = 500]) of the hemoglobin model (20)
14
Parameters
Units
Estimates
Vmax
IU/day
6,890a
Km
IU
1,042a
13
Model parameters 12
11
10
9 0
50
100
150
200
250
300
350
400
450
500
Time t (day)
Fig. 4 Fit of the hemoglobin data with the approximation (trapezoidal rule with n = 500) of the solution (22) of the distributed lifespan model with precursor (20)
HðtÞ ¼ KH MðtÞ;
ð21Þ
where M(t) describes the RBC and H(t) the Hgb concentration. The value KH = 29.5 g/dL is the average amount of Hgb per RBC. The model parameters are h ¼ ðVmax ; Km ; Hen ; Smax ; SC50 ; TN ; bÞ: where Hen is the endogenous level of Hgb due to Een, see [8]. In [8] a three dimensional ordinary differential equation system for (20) based on (18) was constructed with constraints between initial values and model parameter. This reformulation was applied in [8] to fit the Hgb data. Unfortunately, in their reformulation of (20) the initial value is only formulated as a constraint between the initial values of the additional ODEs and not explicitly given. In DLSMs with known inflow for t B 0, according to (7)–(9) there is no freedom to choose 0 M(0). Therefore, due to (9) we have M(0) = kin TM. In this work we apply the solution representation (13) of the precursor model (7), (8) to (20) and obtain Z1 MðtÞ ¼ ð1 LM ðsÞÞkin ðt TN sÞds; ð22Þ 0
where LM(s) is given by (19). The solution (22) is approximated by the trapezoidal rule with n = 500 and a = 0.01. In addition, we apply the PK averaging technique presented in [8] because c(t) has an impulse-like behavior (186 doses over 460 days with rapid clearance). Therefore, Nichols et al. suggested to average the concentration in certain time intervals by Smax h Km Een kin ðdi Þ ¼ Km SC50 Ti Vmax SC50 þ Een SC50 þ Een þ di ln SC50 þ Een i Een Smax Een þ 1 ; ð23Þ di þ SC50 þ Een SC50 þ Een
Hen (CV%)
g/dL
10.26 (1.97)
Smax
cell/day
0.207E-1a
SC50 (CV%)
IU
86.71 (83.1)
TN (CV%) b (CV%)
day
6.98 (10.0) 0.605E-1 (10.0)
TM (CV%)
day
33.1 (10.0)
r2M
day2
547 (20.0)
Secondary parameters (CV%)
Fitting information
a
SofS
68.84
Iter
129
fixed to the values presented in [8]
where Ti is a time interval and di the dose at the time point ti. However, in the following some discussion is necessary. In [8] it is not obvious how Ti was set in (23). Here we choose the interval length Ti for a given dose di as the difference between the time point of the next dose di?1 and the time point of di. Also the PK parameters (Vmax,Km) (no PK data available) are estimated just from the Hgb PD data. This could lead to an over-parameterization of the model. Therefore, we fix Vmax, Km and also Smax to the values presented in the original work. In Fig. 4 the fit of the approximation of (22) with the trapezoidal rule (n = 500) is presented and in Table 4 the parameter estimates of (SC50, Hen, TN, b) are shown from the ADAPT implementation. The fit and the estimates do not exactly coincide due to the differences with the setting of the initial value M(0). Unfortunately, in [8] no derivation of their reformulation of (20) is presented and no information about the optimization process is given.
Discussion Previously published population models with a distributed lifespan (see [1]) based on the convolution operator were mathematically observed. The appearance of the convolution in the rate of change formulation makes an implementation extremely difficult in standard PKPD software. In the presented work we calculated the integral solution representation of DLSMs with and without a precursor population for arbitrary distributions and approximated these representations e.g. by the trapezoidal rule. Therefore, these models could now be simply implemented
123
646
J Pharmacokinet Pharmacodyn (2013) 40:639–650
in every PKPD software. In addition, the presented results are formulated in the general case of a non-constant past for the production kin(t). Further we investigated computational properties of different model formulations. We demonstrated that approximating the calculated solution representation instead of solving the rate of change formulation or the approximation presented in [1], results in a tremendous computational time advantage. Moreover, it was observed that our approach avoids numerical optimization problems. Finally, we applied our method to a published Hgb model and data from [8]. For the sake of completeness we mention that if the production kin(t) is not explicitly given, the advantage of our method will decrease because an additional ODE has to be simultaneously solved. This leads to a relevant increase in computational time. To overcome this issue other approximation techniques are necessary. For the special case of lifespan models with a point distribution it is already proven that an approximation with transit compartments exists, see [15, 16]. However, for lifespan models with an arbitrary continuous distribution an approximation by transit compartments is topic of our next research.
Then lim NðtÞ ¼ 0
t!1
, Nð0Þ ¼
Z1
lN ðxÞ
Z1
0
kf in ðs xÞdsdx:
0
Using (24) Nð0Þ ¼
¼
Z1
lN ðxÞ
Zx
0
0
Z1
Z0
lN ðxÞ
kin ðs xÞdsdx
kin ðsÞdsdx:
x
0
Appendix 2 0 We consider a constant past kin(t) = kin for t B 0 and obtain for the initial values
Nð0Þ ¼
Acknowledgments The present project is supported by the National Research Fund, Luxembourg, and cofunded under the Marie Curie Actions of the European Commission (FP7-COFUND).
Z1
Z0
lN ðxÞ
0 kin dsdx
x
0 0 ¼ kin
Z1
0 lN ðxÞxdx ¼ kin TN ;
0
Appendices
and
Appendix 1 Mð0Þ ¼ To estimate the amount of individuals at the initial population state N(0) for the DLSM (5), (6) we define the production term kin ðtÞ t 0; f kin ðtÞ ¼ ð24Þ 0 t [ 0: Then the appropriate initial value N(0) is determined by the claim limt!1 NðtÞ ¼ 0: Substituting (24) into (5) results in d NðtÞ ¼ kf in ðtÞ dt
Z1
kf in ðt xÞlN ðxÞdx:
ð25Þ
¼
¼
Z1
lM ðyÞ
Z0
0
y
Z1
Z0
lM ðyÞ
0
y
Z1
Z0
lM ðyÞ
0 kin lN ðsÞdsdy
0 kin
Z1
lN ðxÞdxdsdy
0 0 0 kin dsdy ¼ kin TM :
y
0
Appendix 3
0
To estimate the amount M(0) of the follow-up population for the DPLSM (7), (8) a similar approach as in Appendix 1 is applied. We define the production
The general solution of (25) reads NðtÞ ¼ Nð0Þ þ
¼ Nð0Þ
Zt 0 Z1 0
123
kf in ðsÞds
Z t Z1 0
lN ðxÞ
Zt 0
kf in ðs xÞlN ðxÞdxds
0
kf in ðs xÞdsdx:
kf in ðtÞ ¼
ðkin lÞðtÞ 0
t 0; t [ 0;
ð26Þ
and use the condition limt!1 MðtÞ ¼ 0 to determine M(0). Substituting (26) into (7) results in
J Pharmacokinet Pharmacodyn (2013) 40:639–650
MðtÞ ¼ Mð0Þ þ
¼ Mð0Þ
Zt
kf in ðsÞds
0 Z1
Z t Z1 0
lM ðyÞ
0
Zt
647
Appendix 5 kf in ðs yÞlM ðyÞdyds Generally, for an arbitrary continuous function f, the trapezoidal rule reads
0
kf in ðs yÞdsdy;
Zsend
0
f ðsÞds
s0
and therefore, Mð0Þ ¼
¼
Z1
lM ðyÞ
Z1
0
0
Z1
Z0
lM ðyÞ
with the equidistant partition sj = s0 ? (j - 1)h with h ¼ for j = 1,…,m ? 1 of the lifespan interval [s0, send].
kf in ðs yÞdsdy ð27Þ ðkin lN ÞðsÞdsdy:
Appendix 4
Appendix 6 Differentiation of (13) with respect to time t gives Z1 Z1 d 0 MðtÞ ¼ ð1 LM ðxÞÞ kin ðt x sÞlN ðsÞdsdx dt 0
Differentiation of (10) with respect to time t gives d NðtÞ ¼ dt
Z1
¼
0
Z1
lN ðsÞ
0
ð1
0 LN ðsÞÞkin ðt
ð28Þ
send s0 m
y
0
m hX ð f ðsiþ1 Þ þ f ðsi ÞÞ; 2 i¼1
Z1
0
ð1 LM ðxÞÞkin ðt x sÞdxds:
0
ð29Þ
sÞds
0
¼ ½ð1 LN ðsÞÞkin ðt sÞs¼1 s¼0 Z1 ðlN ðsÞÞðkin ðt sÞÞds
The inner integral in (29) could be written as Z1
0
ð1 LM ðxÞÞkin ðt x sÞdx
0 0
¼ kin ðtÞ
Z1
lN ðsÞkin ðt sÞds
0
With (10) we obtain for the initial value Nð0Þ ¼
2 ¼ 4ð1 LN ðsÞÞ
Zs
Z1
ðlN ðsÞÞ
0
¼
0
kin ðt x sÞlM ðxÞdx:
ð31Þ
lN ðsÞ
Zs
Substitution of (31) into (29) gives 3s¼1
kin ðxÞdx5
0
Z1
¼ kin ðt sÞ
Z1 0
ð1 LN ðsÞÞkin ðsÞds
0
ð30Þ
0
¼ kin ðtÞ ðkin lN ÞðtÞ:
Z1
¼ ½ð1 LM ðxÞÞkin ðt x sÞx¼1 x¼0 Z1 kin ðt x sÞlM ðxÞdx
d MðtÞ ¼ dt
Z1 0
s¼0
kin ðxÞdxds
lN ðsÞkin ðt sÞds
Z1 0
lN ðsÞ
Z1
kin ðt x sÞlM ðxÞdxds
0
0
Z0 s
¼ ðkin lN ÞðtÞ kin ðxÞdxds:
Z1
lN ðsÞððkin lM Þðt sÞÞds
0
¼ ðkin lN ÞðtÞ ðkin lN lM ÞðtÞ:
123
648
J Pharmacokinet Pharmacodyn (2013) 40:639–650
For the initial value we obtain Mð0Þ ¼
Z1
ð1 LM ðxÞÞ
Z1
0
Z x Z1
¼ 4ð1 LM ðxÞÞ
0
þ
Z1
lM ðxÞ
0
¼
¼
kin ðx sÞlN ðsÞdsdx
0
2
Z1
Z x Z1 0
lM ðxÞ
Zx
kin ðy sÞlN ðsÞdsdy5
0
x¼0
kin ðy sÞlN ðsÞdsdydx
ðkin lN ÞðyÞdydx
0
0
Z0
0
3x¼1
0
Z1
lM ðxÞ
CDF L(x) of the WBD:
ðkin lN ÞðyÞdydx:
x
Appendix 7 In the following we present the implementation of the trapezoidal rule for (10). Only additional functions and modified ADAPT subroutines are presented. For the complete code see the supplemental material. Logarithm of the gamma function CðxÞ based on [14]:
123
Multiple dosing function for c(t):
J Pharmacokinet Pharmacodyn (2013) 40:639–650
Production kin(t):
649
The ADAPT subroutine OUTPUT calls the trapezoidal rule for the different dosing groups.
Integrand in (10):
Trapezoidal rule:
123
650
References 1. Krzyzanski W, Woo S, Jusko WJ (2006) Pharmacodynamic models for agents that alter production of natural cells with various distributions of lifespans. J Pharmacokinet Pharmacodyn 33(2):125–166
123
J Pharmacokinet Pharmacodyn (2013) 40:639–650 2. Freise K, Widness J, Schmidt R, Veng-Pedersen P (2008) Modeling time variant distributions of cellular lifespans: increases in circulating reticulocyte lifespans following double phlebotomies in sheep. J Pharmacokinet Pharmacodyn 35(3):285–323 3. Ramakrishnan R, Cheung WK, Wacholtz MC, Minton N, Jusko WJ (2004) Pharmacokinetic and pharmacodynamic modeling of recombinant human erythropoietin after single and multiple doses in healthy volunteers. J Clin Pharmacol 44:991–1002 4. Ait-Oudhia S, Scherrmann JM, Krzyzanski W (2010) Simultaneous pharmacokinetics/pharmacodynamics modeling of recombinant human erythropoietin upon multiple intravenous dosing in rats. J Pharmacol Exp Ther 334:897–910 5. D’Argenio DZ, Schumitzky A, Wang X (2009) ADAPT 5 user’s guide: pharmacokinetic/pharmacodynamic systems analysis software. Biomedical Simulations Resource, Los Angeles 6. Beal S, Sheiner LB, Boeckmann A, Bauer RJ (2009) NONMEM user’s guides. Icon Development Solutions, Ellicott City 7. MATLAB Release (2012) The MathWorks, Inc., Natick 8. Nichols B, Shrestha RP, Horowitz J, Hollot CV, Germain MJ, Gaweda AE, Chait Y (2011) Simplification of an erythropoiesis model for design of anemia management protocols in end stage renal disease. In: 33rd Annual international conference of the IEEE EMBS Boston, Massachusetts, USA 9. Atkinson K (1989) An introduction to numerical analysis, 2nd edn. Wiley, New York 10. Thomas GB, Finney RL (1996) Calculus and analytic geometry, 9th edn. Addison-Wesley, Reading 11. Marshall AW, Olkin I (2007) Life distributions. Springer, New York 12. Bebbinton M, Lai CD, Zitikis R (2007) Modeling human mortality using mixtures of bathtub shaped failure distributions. J Theor Biol 245(3):528–538 13. Korell J, Coulter CV, Duffull SB (2011) A statistical model for red blood cell survival. J Theor Biol 268(1):39–49 14. Press WH, Teukolsky SA, Vetterling WT, Flannery BP (1992) Numerical recipes in Fortran 77, the art of scientific computing, 2nd edn. Cambridge University Press, New York 15. Krzyzanski W (2011) Interpretation of transit compartments pharmacodynamic models as lifespan based indirect response models. J Pharmacokinet Pharmacodyn 38(2):179–204 16. Koch G, Schropp J (2012) General relationship between transit compartments and lifespan models. J Pharmacokinet Pharmacodyn 39(4):343–355
