- Email: [email protected]
JANA: A new iterative polyexponential curve stripping program
Computer Methods and Programs in Biomedicine 20 (1985) 269-275
269
Elsevier CPB 00710
JANA: A new iterative polyexponential curve stripping program Adrian Dunne Department of Pharmacology, University College Dublin, Belfield, Dublin 4, Ireland
A new iterative polyexponential curve stripping technique for the provision of initial pharmacokinetic parameter estimates has been developed. Such estimates are required for nonlinear least-squares curve fitting. In contrast to conventional curve stripping, this new technique does not make any assumption about the relative magnitudes of the exponential rate constants. Hence, the parameter estimates which it provides are free of the bias which may arise in conventional curve stripping. A BASIC program called JANA has been developed to implement the new curve stripping procedure. Curve stripping
Curve fitting
Pharmacokinetic parameter estimation
1. Introduction The pharmacokinetic characteristics of a drug are frequently described by polyexponential models of the form: m
C= Y'~ Ag.eXp(-Bj.t)
(1)
j=l
where C is the plasma drug concentration at time t following administration, m is the number of exponential terms and Aj and Bj are the coefficients and exponents [1]. Estimates of the parameters (Aj and Bj) can be obtained from experimental data using a curve fitting or parameter estimation technique. The most widely used procedure is the nonlinear weighted least-squares method [2,3] which uses a computerised search to locate the minimum on the so-called 'weighted residual sum of squares surface'. A number of iterative search algorithms are available [4] and some of these have been used as the basis for computer programs for nonlinear weighted least-squares estimation of pharmacokinetic parameters [5-10]. All of these algorithms require initial estimates of the parameters in order to begin the search and successful
termination of the algorithm may depend on how good the initial parameter estimates are. Pedersen [8] has said 'regardless of the choice of algorithm, the question of which starting values the parameters should be given still seems to be the greatest practical problem the user faces in nonlinear estimation'. Several graphical/numerical methods have been employed to generate initial pharmacokinetic parameter estimates [11-17]. The most popular of these has been curve stripping, which is also known as back projection, the method of residuals, feathering or peeling. This method assumes that the exponential terms have disparate rate constants and has been automated for execution by digital computers [18], microcomputers [19,20] and hand held calculators [21]. However, attention has been drawn to the fact that curve stripping produces biased parameter estimates if the assumption regarding the relative magnitudes of the rate constants is invalid [3,8,22]. JANA is a program based on a modified curve stripping technique which does not make any assumption about the relative magnitudes of the rate constants. It is written in BASIC and may be implemented on a microcomputer.
270
2. Theory Consider the situation in which a pharmacokinetic study yields n data pairs ( C i, t i i = 1 . . . . . n ) where C'~ represents the measured plasma drug concentration at time t i following administration. Initial estimates of the parameters ( A j, Bj) of eq. (1) are required for nonlinear least-squares curve fitting.
If it is now assumed that there exists a time T2 such that for t > T 2 the first m - 2 exponential terms have decayed, i.e. m--2
Y~ A j " e x p ( - B j , t) = 0
(8)
then R, = A , , _ , " e x p ( - B m _ t "t)
2.1. Conventional curve stripping
t > T2
j=l
t > T2.
(9)
Taking natural logarithms this yields Let the exponents be ordered such that (2)
B 1 > B 2 > B 3 > ..... > Bm
and assume that there exists a time T~ such that for t > T 1 the first m - 1 exponential terms have decayed, i.e.: rrt--I
E
Aj" e x p ( - B j , t ) = 0
t > T 1,
(3)
j=l
then C=A,..exp(-B.,.t)
(4)
t > T 1.
In R t = l n A , , _ l - B i n _
,'t
A
m
-
Bm • t
-
RE,i=Rl,i-,4m_~'exp(-Bm_l"ti)
which predicts that for t > Tt a semilogarithmic plot of C versus t should be linear with slope - B , , and intercept In A m. Hence, if the observations (C i, t/) are plotted semilogarithmically the terminal points may be used to estimate A,, and B m by linear least-squares regression. These estimates are denoted by Am and Bin- The values of the last exponential term at each of the earlier time points may be estimated, using "4m and Bm, and subtracted from the observed plasma concentrations yielding the first set of residuals Rz,i=Ci-Am'eXp(-hm'ti)
t i < T ,.
(6)
The model for these residuals consists of the first m - 1 exponential terms, i.e.: m--1
R , = Y'. A j . e x p ( - B j . t ) . j=l
ti
(11) (5)
t > 7"1,
(10)
which predicts that for t > T2 a semilogarithmic plot of R 1 versus t should be linear with slope B,,,_ 1 and intercept In A,,_ p Consequently, the terminal points of a semilogarithmic plot of the residuals (R], i) versus time (ti) may be used to produce estimates A.,_ ] and Bin-1 of A . , _ ] and B.,_~ by linear least-squares regression. Proceeding in this manner a second set of residuals
Taking natural logarithms, this yields In C = In
t > T 2,
(7)
is computed and the terminal points used to estimate Am_ 2 and Bin_ 2. This process is repeated until all the parameters have been estimated. The curve stripping technique is illustrated graphically in Fig. 1. ,~ too' _a
O
t~ 0.10
Q < ~ o.ol TIME
Fig. 1. Three exponential terms being stripped from data collected following bolus intravenous drug administration, x , Plasma drug concentrations; e, first set of residuals; O, second set of residuals.
271
If the data were collected following extravascular drug administration then a lag time must be considered. The presence of a lag time is indicated if
Aj <
(12)
O.
j~l
An estimate of the lag time ( L ) is found by trial and error solution of ~.~ Aj. e x p ( - B j . L ) = 0
(13)
j=l
and the coefficients adjusted accordingly. The above procedure is based on the assumption that there exists T, (k = 1. . . . . m - 1) such that for t > T k the first m - k exponential terms have decayed, i.e.: m-k
Y'~ A j . e x p ( - B j . t ) = O
t > T k.
(14)
j=l
This assumption may be approximately satisfied if: B1 >> B2
>>
B 3 >>
.....
>> B m .
(15)
If this condition is not met the curve stripping procedure will yield biased estimates of the parameters. 2.2. Modified curve stripping
Assumption (14) would be unnecessary if for t > Tk m-k
the term
~ A j - e x p ( - B j . t) could be evaluated j=l
at the sample times (ti) and subtracted from the corresponding Rk_l. i (or Ci when k = 1) to produce 'corrected' residuals m-k R *, - , , , =
R*-l,i-
E Aj e x p ( - B j , ti).
(16)
j=l
or 'corrected' plasma concentrations m-1
Ci* = C i -
Y'~ A j " e x p ( - B / , t,). j~l
(17)
Curve stripping using 'corrected' plasma concentrations or residuals would produce unbiased parameter estimates. However, the corrections applied in eqs. (16) and (17) are impossible in practice since they require a knowledge of the parameters Aj and Bj, which are to be estimated. Instead, initial estimates of the parameters (obtained by making assumption (14)) are used in place of the parameters themselves in expressions (16) and (17) and a second set of parameter estimates calculated using the 'corrected' data. Since it can be shown that these latter parameter estimates are better than the initial ones, the method is applied iteratively as follows: (1) Obtain parameter estimates by conventional curve stripping. (2) Apply 'corrections' (16) and (17) using these estimates in place of the parameters. (3) Using 'corrected' plasma concentrations and residuals obtain parameter estimates by curve stripping. (4) Repeat steps (2) and (3) until the fit of the model to the data improves by less than a specified fraction.
3. Program description The overall structure of JANA is outlined in the flowchart in Fig. 2. It begins with one or two exponential terms (depending on the route of drug administration) and increases in steps of one up to the maximum number of exponential terms specified by the user. Presently the program can fit a maximum of four exponential terms; extending it to fit five or more terms should be straightforward. The program considers every possible partition of the data among the exponential terms, with a minimum of two points in each exponential phase. That partition which has the smallest weighted residual sum of squares is selected. For data collected following extravascular drug administration the possibility of a lag time must be considered. The procedure for computing a lag time is outlined in Fig. 3. If the lag time estimate is less than 25% of the first nonzero sample time then the lag time is included provided i t results in a better fit of the data than forcing the fitted curve
272
SPECIFICATIONS l
I COMPUTE LAGTIME BY I FIT THROUGH I
[ TRIAL AND ERROR
J
I
I
ORIGIN
NO t YES
I
.,-,
I
JCALCULATE W.RSS.ANDI STORE COEFFICIENTS [ I
JFIT THROUGH ORIGIN
USING UNADJUSTED COEFFICIENTS
NO+ i
J
CALCULATE W.R~.S.
YES . COMPUTE LAG TIME RETRIEVE ADJUSTED J COEFFICIENTS FROM STORAGE
J
J
PRINT RESULTS LT= 1
-I
CONTINUE
J I
Fig. 3. Flowchart for lag time estimation. SP is the sum of the exponential coefficients, LT is a flag which indicates whether or not a lag time is required, TL is the lag time estimate, TT is the first non-zero sample time and W.R.S.S. is the weighted residual sum of squares.
Fig. 2. Flowchart showing the overall structure of the program J A N A . N P is the number of parameters in the model, i.e. twice the number of exponential terms. Extravascular drug administration is indicated by RT = 1. MX is the maximum number of exponential terms of interest.
through the origin. When a lag time is not indicated the fitted curve is forced through the orgin. The program offers a menu of five alternative weighting schemes. The weights may be input with the data or generated by the program using W~ = 1 / G q
(18)
273
CLii¢¢E
I i [Lt
z
I}i
COI.iVE b ' G E i q ( L
U t < i I EF( [ I } N
F i Lh
Ei:'-~:i;/ i l ~
I W]
L]~:
F:'t~B:,qME I El<
i ~ L : R A I iL!Hb
A
=
ot~.o.Zli~x.i4
L;
=
06.03410o4
E: =
i.A U
(I ~
i: ,,f:'[_l~}6i,~t J ~,i ':
.....
F i i i E:i)
_,2
=
1.V97275/5 .bi229,q*~41
NUMBF:i;: NUME~EP.
I.]F
F'LI i b,I I c:.'~ =
3
OF
POLNIS
6
I 1HE
=
] .Z.::ti. ~ 0
=
15. 590441,~J
=
.9U?30~EiSii~
L:-L)LIS
IL:-(]AL(I
[[}~IS'IgALC,
i_/
i )
2
. 1662
,~; 4 5; 6 2 8 u,
. :2; &:~;:i!!; ,, 5 I Z # ,5 8
15. / J~:i:.!;. i 131 . :d 3t . ~ 2 8 , :d
0 14, 3090~349 .:' :2;. ~., J i ), 2 9 Z 9 . : ;"Y.21;6~J :; 'L";Z!. 7 4 0 J ~"4:3 Z 7 . M ~ 5 : S q i ~.: 12. (~4q5] O, 5.54.;~ 59858 :,.:. 4 : ~ : : 9 8 0 5 0 5
li J. S'?~;:,'1511 .... . t5 J (.)02[.~'~* : 3 i ~ 8 0 ~ : 5 1 :~;I ~q :i 4 0 1 2 4 3 3 . "~ 9 4 0 5 ~ . ~ 8 5 8 . i5548%50:2!; .052(:,,)i 4 175 . ()Z:'~OO504,b5
J
=
f{FL!Uih:EI)
liME
E:UI~F
I:'EilNi
IL
L S 1 I I"l~q ] E!5
A
iqU
the maximum number of iterations required and the convergence criterion. The parameter estimates are subject to certain constraints: only negative exponentials are considered and there is no restriction on the ratio of the exponents. Furthermore, when the data correspond with intravenous drug administration the coefficients must all be positive. Data collected after intramuscular administration of spectinomycin have been used for testing curve stripping programs [18,19]. These data were analysed using J A N A and the output is shown in Fig. 4. The program is written in Applesoft BASIC for implementation on an Apple lie 64K microcomputer with an 80 column card and an Epson RX-80 printer. The disk operating system is DOS 3.3 and the program occupies 11.3 Kbytes of memory.
[Iqb
(]1 I } q H M Y ( Z i N
S['t
i)(-, ! ¢:, t: [ U R E I )
NUFIUEF;:
b I t<] t t -
IZ.8 5.0
:2. 4
Wh]:GHf 1 1 J ] J J ]
4. Results
1.
The conventional curve stripping program CSTRIP [18] and J A N A were compared using simulated pharmacokinetic data. Noise-free data were generated using a triexponential model corresponding to intravenous drug administration, i.e.:
J
Fig. 4. Sample output from JANA.
where q may take one of the values 0, 0.5, 1 or 2. The weights are normalised by the program so that their sum is equal to the total nuhaber of observations. The user is asked to specify the maximum number of exponential terms of interest, the route of drug administration (intravenous or extravascular),
y
=A 1 • exp(-B,-
t) + A 2 • e x p ( - B 2 • t)
+ A 3 • e x p ( - B 3 - t)
(19)
TABL E 1 Parameter estimates produced by CSTRIP and JANA using simulated data a Experiment
Parameter A1
B1
A2
B2
A3
B3
1
True CSTRIP JANA
33.33 25.85 33.34
0.10 0.137 0.100
33.33 38.59 33.33
0.02 0.026 0.020
33.33 36.50 33.33
0.004 0.0042 0.0040
2
True CSTRIP JANA
33.33 9.12 33.33
0.06 0.123 0.060
33.33 49.76 33.33
0.02 0.034 0.020
33.33 41.51 33.34
0.007 0.0075 0.0070
a Noise free data were generated from a triexponential (intravenous drug administration) model with the parameter values labelled ' true'.
274 TABLE 2 P a r a m e t e r estimates p r o d u c e d b y C S T R I P a n d J A N A using s i m u l a t e d d a t a a Experiment
Parameter A1
B1
A2
B2
L
3
True CSTRIP JANA
- 50.00 - 48.95 - 50.10
0.10 0.105 0.100
50.00 48.88 50.02
0.02 0.020 0.020
5.0 5.14 4.98
4
True CSTRIP JANA
- 50.00 - 27.58 - 49.98
0.04 0.058 0.040
50.00 27.56 49.96
0.02 0.016 0.020
5.0 5.52 4.98
a N o i s e free d a t a were g e n e r a t e d from a b i e x p o n e n t i a l ( e x t r a v a s c u l a r d r u g a d m i n i s t r a t i o n ) m o d e l with the p a r a m e t e r values labelled • true'.
with A1 = A2 = A3 = 33.33, B1 = 0.10, B2 = 0.02 and B 3 = 0.004. There were twenty data values at t = 0, 1, 2, 3, 5, 7, 10, 15, 20, 25, 30, 40, 50, 70, 100, 150, 200, 250, 350 and 500 units. Three exponential terms were stripped from the data using CSTRIP and JANA and the results are tabulated in Table 1 (experiment 1). A second set of data was generated at the same sample times using the same model with the ratios of the exponential rate constants decreased. This data was analysed by CSTRIP and JANA and the results are also summarised in Table 1 (experiment 2). Two further noise-free data sets were generated using the biexponential model corresponding to extravascular drug administration with a lag phase, i.e.:
stripping technique does not depend on condition (15) and consequently its parameter estimates were not subject to the same degree of bias. JANA begins with a conventional curve stripping of the data and then proceeds to improve on the parameter estimates in the manner described above. Therefore, at worst it gives parameter estimates identical with those produced by conventional curve stripping. At best JANA provides parameter estimates which are as good as those arrived at by nonlinear least-squares regression. In general it provides better initial estimates than conventional curve stripping and thus increases the likelihood of a successful nonlinear leastsquares fit, with fewer iterations being required.
y = A 1. e x p ( - B , ( t - L ) )
6. Availability
+A 2 • e x p ( - B 2 ( t - L ) ) .
(20)
There were fifteen samples at t = 5, 6, 7, 8, 9, 10, 12, 15, 20, 25, 35, 55, 80, 105 and 155 units. These data sets were analysed by CSTRIP and JANA and the results are tabulated in Table 2.
5. Discussion Tables 1 and 2 show that as condition (15) above became less tenable the conventional curve stripping procedure produced increasingly biased parameter estimates. However, the modified curve
Details of the program's availability can be obtained from the author.
Acknowledgements The author wishes to acknowledge the support of Elan Corporation p.l.c, and the Medical Research Council of Ireland.
References [1] M. G i b a l d i a n d D. Perrier, P h a r m a c o k i n e t i c s (Marcel Dekker, N e w York, 1982).
275 [2] B.E. Rodda, Analysis of sets of estimates from pharmacokinetic studies, in: Kinetic Data Analysis, ed. L. Endrenyi, pp. 285-298 (Plenum Press, New York, 1981). [3] "C.C. Peck and B.B. Barrett, Nonlinear least-squares regression programs for microcomputers, J. Pharmacokin. Biopharm. 7 (1979) 537-541. [4] Y. Bard, Nonlinear Parameter Estimation (Academic Press, New York, 1974). [5] M. Pfeffer, COMPT, a time sharing program for nonlinear regression analysis of compartmental models of drug distribution, J. Pharmacokin. Biopharm. 1 (1973) 137-163. [6] C.M. Metzler, G.L. Elfring and A.J. McEwen, A package of computer programs for pharmacokinetic modeling, Biometrics 30 (1974) 562. [7] J.G. Wagner, Fundamentals of Clinical Pharmacokinetics (Drug Intelligence, Illinois, 1975). [8] P.V. Pedersen, Curve fitting and modelling in pharmacokinetics and some practical experiences with NONLIN and a new program FUNFIT, J. Pharmacokin. Biopharm, 5 (1977) 513-531. [9] R. Gomeni and C. Gomeni, AUTOMOD: A polyalgorithm for an integrated analysis of linear pharmacokinetic models, Comput. Biol. Med. 9 (1979) 39-48. [10] A. Messori, G.D. Cori and E. Tendi, Iterative least-squares fitting programs in pharmacokinetics for a programmable hand held calculator, Am. J. Hosp. Pharm. 40 (1983) 1673-1684. [11] D.G. Gardner, J.C. Gardner and G. Laush, Method for the analysis of multicomponent exponential decay curves, J. Chem. Phys. 31 (1959) 978-986. [12] M.R. Smith and S.T. Nichols, Improved resolution in the analysis of multicomponent exponential signals, Nuclear lnstrum. Meth. 205 (1983) 479-483.
[13] R.G. Cornell, A method for fitting linear combinations of exponentials, Biometrics 18 (1962) 104-113. [14] D.H. Parsons, Biological problems involving sums of exponential functions of time: an improved method of calculation, Math. Biosci. 9 (1970) 37-47. [15] S.D. Foss, A method of exponential curve fitting by numerical integration, Biometrics 26 (1970) 815-821. [16] K.C. Yeh and K.C. Kwan, Kinetic parameter estimation by numerical algorithms and multiple linear regression: application to pharmacokinetics, J. Pharm. Sci. 68 (1979) 1120-1124. [17] J.R. Koup, Direct linear plotting method for estimation of pharmacokinetics parameters, J. Pharm. Sci. 70 (1981) 1093-1094. [18] A.J. Sedman and J.G. Wagner, CSTRIP, a FORTRAN IV computer program for obtaining polyexponential parameter estimates. J. Pharm. Sci. 65 (1976) 1006-1010. [19] R.D. Brown and J.E. Manno, ESTRIP: a BASIC computer program for obtaining initial polyexponential parameter estimates, J. Pharm. Sci. 67 (1978) 1687-1691. [20] J.G. Leferink and R.A.A. Maes, STRIPACT, an interactive curve fit programme for pharmacokinetic analyses. Arzneim. Forsch. 29 (1979) 1894-1898. [21] P.P. leBlanc and J. Dumas, Calcul des valeurs initiales des parametres pharmacocinetiques a raide d'un calculateur programmable, Therapie 38 (1983) 21-26. [22] A. Dunne and L. Wilson, Multiexponential parameter estimation: elimination of bias, Br. J. Pharmacol. 80 (1983) 714P.
269
Elsevier CPB 00710
JANA: A new iterative polyexponential curve stripping program Adrian Dunne Department of Pharmacology, University College Dublin, Belfield, Dublin 4, Ireland
A new iterative polyexponential curve stripping technique for the provision of initial pharmacokinetic parameter estimates has been developed. Such estimates are required for nonlinear least-squares curve fitting. In contrast to conventional curve stripping, this new technique does not make any assumption about the relative magnitudes of the exponential rate constants. Hence, the parameter estimates which it provides are free of the bias which may arise in conventional curve stripping. A BASIC program called JANA has been developed to implement the new curve stripping procedure. Curve stripping
Curve fitting
Pharmacokinetic parameter estimation
1. Introduction The pharmacokinetic characteristics of a drug are frequently described by polyexponential models of the form: m
C= Y'~ Ag.eXp(-Bj.t)
(1)
j=l
where C is the plasma drug concentration at time t following administration, m is the number of exponential terms and Aj and Bj are the coefficients and exponents [1]. Estimates of the parameters (Aj and Bj) can be obtained from experimental data using a curve fitting or parameter estimation technique. The most widely used procedure is the nonlinear weighted least-squares method [2,3] which uses a computerised search to locate the minimum on the so-called 'weighted residual sum of squares surface'. A number of iterative search algorithms are available [4] and some of these have been used as the basis for computer programs for nonlinear weighted least-squares estimation of pharmacokinetic parameters [5-10]. All of these algorithms require initial estimates of the parameters in order to begin the search and successful
termination of the algorithm may depend on how good the initial parameter estimates are. Pedersen [8] has said 'regardless of the choice of algorithm, the question of which starting values the parameters should be given still seems to be the greatest practical problem the user faces in nonlinear estimation'. Several graphical/numerical methods have been employed to generate initial pharmacokinetic parameter estimates [11-17]. The most popular of these has been curve stripping, which is also known as back projection, the method of residuals, feathering or peeling. This method assumes that the exponential terms have disparate rate constants and has been automated for execution by digital computers [18], microcomputers [19,20] and hand held calculators [21]. However, attention has been drawn to the fact that curve stripping produces biased parameter estimates if the assumption regarding the relative magnitudes of the rate constants is invalid [3,8,22]. JANA is a program based on a modified curve stripping technique which does not make any assumption about the relative magnitudes of the rate constants. It is written in BASIC and may be implemented on a microcomputer.
270
2. Theory Consider the situation in which a pharmacokinetic study yields n data pairs ( C i, t i i = 1 . . . . . n ) where C'~ represents the measured plasma drug concentration at time t i following administration. Initial estimates of the parameters ( A j, Bj) of eq. (1) are required for nonlinear least-squares curve fitting.
If it is now assumed that there exists a time T2 such that for t > T 2 the first m - 2 exponential terms have decayed, i.e. m--2
Y~ A j " e x p ( - B j , t) = 0
(8)
then R, = A , , _ , " e x p ( - B m _ t "t)
2.1. Conventional curve stripping
t > T2
j=l
t > T2.
(9)
Taking natural logarithms this yields Let the exponents be ordered such that (2)
B 1 > B 2 > B 3 > ..... > Bm
and assume that there exists a time T~ such that for t > T 1 the first m - 1 exponential terms have decayed, i.e.: rrt--I
E
Aj" e x p ( - B j , t ) = 0
t > T 1,
(3)
j=l
then C=A,..exp(-B.,.t)
(4)
t > T 1.
In R t = l n A , , _ l - B i n _
,'t
A
m
-
Bm • t
-
RE,i=Rl,i-,4m_~'exp(-Bm_l"ti)
which predicts that for t > Tt a semilogarithmic plot of C versus t should be linear with slope - B , , and intercept In A m. Hence, if the observations (C i, t/) are plotted semilogarithmically the terminal points may be used to estimate A,, and B m by linear least-squares regression. These estimates are denoted by Am and Bin- The values of the last exponential term at each of the earlier time points may be estimated, using "4m and Bm, and subtracted from the observed plasma concentrations yielding the first set of residuals Rz,i=Ci-Am'eXp(-hm'ti)
t i < T ,.
(6)
The model for these residuals consists of the first m - 1 exponential terms, i.e.: m--1
R , = Y'. A j . e x p ( - B j . t ) . j=l
ti
(11) (5)
t > 7"1,
(10)
which predicts that for t > T2 a semilogarithmic plot of R 1 versus t should be linear with slope B,,,_ 1 and intercept In A,,_ p Consequently, the terminal points of a semilogarithmic plot of the residuals (R], i) versus time (ti) may be used to produce estimates A.,_ ] and Bin-1 of A . , _ ] and B.,_~ by linear least-squares regression. Proceeding in this manner a second set of residuals
Taking natural logarithms, this yields In C = In
t > T 2,
(7)
is computed and the terminal points used to estimate Am_ 2 and Bin_ 2. This process is repeated until all the parameters have been estimated. The curve stripping technique is illustrated graphically in Fig. 1. ,~ too' _a
O
t~ 0.10
Q < ~ o.ol TIME
Fig. 1. Three exponential terms being stripped from data collected following bolus intravenous drug administration, x , Plasma drug concentrations; e, first set of residuals; O, second set of residuals.
271
If the data were collected following extravascular drug administration then a lag time must be considered. The presence of a lag time is indicated if
Aj <
(12)
O.
j~l
An estimate of the lag time ( L ) is found by trial and error solution of ~.~ Aj. e x p ( - B j . L ) = 0
(13)
j=l
and the coefficients adjusted accordingly. The above procedure is based on the assumption that there exists T, (k = 1. . . . . m - 1) such that for t > T k the first m - k exponential terms have decayed, i.e.: m-k
Y'~ A j . e x p ( - B j . t ) = O
t > T k.
(14)
j=l
This assumption may be approximately satisfied if: B1 >> B2
>>
B 3 >>
.....
>> B m .
(15)
If this condition is not met the curve stripping procedure will yield biased estimates of the parameters. 2.2. Modified curve stripping
Assumption (14) would be unnecessary if for t > Tk m-k
the term
~ A j - e x p ( - B j . t) could be evaluated j=l
at the sample times (ti) and subtracted from the corresponding Rk_l. i (or Ci when k = 1) to produce 'corrected' residuals m-k R *, - , , , =
R*-l,i-
E Aj e x p ( - B j , ti).
(16)
j=l
or 'corrected' plasma concentrations m-1
Ci* = C i -
Y'~ A j " e x p ( - B / , t,). j~l
(17)
Curve stripping using 'corrected' plasma concentrations or residuals would produce unbiased parameter estimates. However, the corrections applied in eqs. (16) and (17) are impossible in practice since they require a knowledge of the parameters Aj and Bj, which are to be estimated. Instead, initial estimates of the parameters (obtained by making assumption (14)) are used in place of the parameters themselves in expressions (16) and (17) and a second set of parameter estimates calculated using the 'corrected' data. Since it can be shown that these latter parameter estimates are better than the initial ones, the method is applied iteratively as follows: (1) Obtain parameter estimates by conventional curve stripping. (2) Apply 'corrections' (16) and (17) using these estimates in place of the parameters. (3) Using 'corrected' plasma concentrations and residuals obtain parameter estimates by curve stripping. (4) Repeat steps (2) and (3) until the fit of the model to the data improves by less than a specified fraction.
3. Program description The overall structure of JANA is outlined in the flowchart in Fig. 2. It begins with one or two exponential terms (depending on the route of drug administration) and increases in steps of one up to the maximum number of exponential terms specified by the user. Presently the program can fit a maximum of four exponential terms; extending it to fit five or more terms should be straightforward. The program considers every possible partition of the data among the exponential terms, with a minimum of two points in each exponential phase. That partition which has the smallest weighted residual sum of squares is selected. For data collected following extravascular drug administration the possibility of a lag time must be considered. The procedure for computing a lag time is outlined in Fig. 3. If the lag time estimate is less than 25% of the first nonzero sample time then the lag time is included provided i t results in a better fit of the data than forcing the fitted curve
272
SPECIFICATIONS l
I COMPUTE LAGTIME BY I FIT THROUGH I
[ TRIAL AND ERROR
J
I
I
ORIGIN
NO t YES
I
.,-,
I
JCALCULATE W.RSS.ANDI STORE COEFFICIENTS [ I
JFIT THROUGH ORIGIN
USING UNADJUSTED COEFFICIENTS
NO+ i
J
CALCULATE W.R~.S.
YES . COMPUTE LAG TIME RETRIEVE ADJUSTED J COEFFICIENTS FROM STORAGE
J
J
PRINT RESULTS LT= 1
-I
CONTINUE
J I
Fig. 3. Flowchart for lag time estimation. SP is the sum of the exponential coefficients, LT is a flag which indicates whether or not a lag time is required, TL is the lag time estimate, TT is the first non-zero sample time and W.R.S.S. is the weighted residual sum of squares.
Fig. 2. Flowchart showing the overall structure of the program J A N A . N P is the number of parameters in the model, i.e. twice the number of exponential terms. Extravascular drug administration is indicated by RT = 1. MX is the maximum number of exponential terms of interest.
through the origin. When a lag time is not indicated the fitted curve is forced through the orgin. The program offers a menu of five alternative weighting schemes. The weights may be input with the data or generated by the program using W~ = 1 / G q
(18)
273
CLii¢¢E
I i [Lt
z
I}i
COI.iVE b ' G E i q ( L
U t < i I EF( [ I } N
F i Lh
Ei:'-~:i;/ i l ~
I W]
L]~:
F:'t~B:,qME I El<
i ~ L : R A I iL!Hb
A
=
ot~.o.Zli~x.i4
L;
=
06.03410o4
E: =
i.A U
(I ~
i: ,,f:'[_l~}6i,~t J ~,i ':
.....
F i i i E:i)
_,2
=
1.V97275/5 .bi229,q*~41
NUMBF:i;: NUME~EP.
I.]F
F'LI i b,I I c:.'~ =
3
OF
POLNIS
6
I 1HE
=
] .Z.::ti. ~ 0
=
15. 590441,~J
=
.9U?30~EiSii~
L:-L)LIS
IL:-(]AL(I
[[}~IS'IgALC,
i_/
i )
2
. 1662
,~; 4 5; 6 2 8 u,
. :2; &:~;:i!!; ,, 5 I Z # ,5 8
15. / J~:i:.!;. i 131 . :d 3t . ~ 2 8 , :d
0 14, 3090~349 .:' :2;. ~., J i ), 2 9 Z 9 . : ;"Y.21;6~J :; 'L";Z!. 7 4 0 J ~"4:3 Z 7 . M ~ 5 : S q i ~.: 12. (~4q5] O, 5.54.;~ 59858 :,.:. 4 : ~ : : 9 8 0 5 0 5
li J. S'?~;:,'1511 .... . t5 J (.)02[.~'~* : 3 i ~ 8 0 ~ : 5 1 :~;I ~q :i 4 0 1 2 4 3 3 . "~ 9 4 0 5 ~ . ~ 8 5 8 . i5548%50:2!; .052(:,,)i 4 175 . ()Z:'~OO504,b5
J
=
f{FL!Uih:EI)
liME
E:UI~F
I:'EilNi
IL
L S 1 I I"l~q ] E!5
A
iqU
the maximum number of iterations required and the convergence criterion. The parameter estimates are subject to certain constraints: only negative exponentials are considered and there is no restriction on the ratio of the exponents. Furthermore, when the data correspond with intravenous drug administration the coefficients must all be positive. Data collected after intramuscular administration of spectinomycin have been used for testing curve stripping programs [18,19]. These data were analysed using J A N A and the output is shown in Fig. 4. The program is written in Applesoft BASIC for implementation on an Apple lie 64K microcomputer with an 80 column card and an Epson RX-80 printer. The disk operating system is DOS 3.3 and the program occupies 11.3 Kbytes of memory.
[Iqb
(]1 I } q H M Y ( Z i N
S['t
i)(-, ! ¢:, t: [ U R E I )
NUFIUEF;:
b I t<] t t -
IZ.8 5.0
:2. 4
Wh]:GHf 1 1 J ] J J ]
4. Results
1.
The conventional curve stripping program CSTRIP [18] and J A N A were compared using simulated pharmacokinetic data. Noise-free data were generated using a triexponential model corresponding to intravenous drug administration, i.e.:
J
Fig. 4. Sample output from JANA.
where q may take one of the values 0, 0.5, 1 or 2. The weights are normalised by the program so that their sum is equal to the total nuhaber of observations. The user is asked to specify the maximum number of exponential terms of interest, the route of drug administration (intravenous or extravascular),
y
=A 1 • exp(-B,-
t) + A 2 • e x p ( - B 2 • t)
+ A 3 • e x p ( - B 3 - t)
(19)
TABL E 1 Parameter estimates produced by CSTRIP and JANA using simulated data a Experiment
Parameter A1
B1
A2
B2
A3
B3
1
True CSTRIP JANA
33.33 25.85 33.34
0.10 0.137 0.100
33.33 38.59 33.33
0.02 0.026 0.020
33.33 36.50 33.33
0.004 0.0042 0.0040
2
True CSTRIP JANA
33.33 9.12 33.33
0.06 0.123 0.060
33.33 49.76 33.33
0.02 0.034 0.020
33.33 41.51 33.34
0.007 0.0075 0.0070
a Noise free data were generated from a triexponential (intravenous drug administration) model with the parameter values labelled ' true'.
274 TABLE 2 P a r a m e t e r estimates p r o d u c e d b y C S T R I P a n d J A N A using s i m u l a t e d d a t a a Experiment
Parameter A1
B1
A2
B2
L
3
True CSTRIP JANA
- 50.00 - 48.95 - 50.10
0.10 0.105 0.100
50.00 48.88 50.02
0.02 0.020 0.020
5.0 5.14 4.98
4
True CSTRIP JANA
- 50.00 - 27.58 - 49.98
0.04 0.058 0.040
50.00 27.56 49.96
0.02 0.016 0.020
5.0 5.52 4.98
a N o i s e free d a t a were g e n e r a t e d from a b i e x p o n e n t i a l ( e x t r a v a s c u l a r d r u g a d m i n i s t r a t i o n ) m o d e l with the p a r a m e t e r values labelled • true'.
with A1 = A2 = A3 = 33.33, B1 = 0.10, B2 = 0.02 and B 3 = 0.004. There were twenty data values at t = 0, 1, 2, 3, 5, 7, 10, 15, 20, 25, 30, 40, 50, 70, 100, 150, 200, 250, 350 and 500 units. Three exponential terms were stripped from the data using CSTRIP and JANA and the results are tabulated in Table 1 (experiment 1). A second set of data was generated at the same sample times using the same model with the ratios of the exponential rate constants decreased. This data was analysed by CSTRIP and JANA and the results are also summarised in Table 1 (experiment 2). Two further noise-free data sets were generated using the biexponential model corresponding to extravascular drug administration with a lag phase, i.e.:
stripping technique does not depend on condition (15) and consequently its parameter estimates were not subject to the same degree of bias. JANA begins with a conventional curve stripping of the data and then proceeds to improve on the parameter estimates in the manner described above. Therefore, at worst it gives parameter estimates identical with those produced by conventional curve stripping. At best JANA provides parameter estimates which are as good as those arrived at by nonlinear least-squares regression. In general it provides better initial estimates than conventional curve stripping and thus increases the likelihood of a successful nonlinear leastsquares fit, with fewer iterations being required.
y = A 1. e x p ( - B , ( t - L ) )
6. Availability
+A 2 • e x p ( - B 2 ( t - L ) ) .
(20)
There were fifteen samples at t = 5, 6, 7, 8, 9, 10, 12, 15, 20, 25, 35, 55, 80, 105 and 155 units. These data sets were analysed by CSTRIP and JANA and the results are tabulated in Table 2.
5. Discussion Tables 1 and 2 show that as condition (15) above became less tenable the conventional curve stripping procedure produced increasingly biased parameter estimates. However, the modified curve
Details of the program's availability can be obtained from the author.
Acknowledgements The author wishes to acknowledge the support of Elan Corporation p.l.c, and the Medical Research Council of Ireland.
References [1] M. G i b a l d i a n d D. Perrier, P h a r m a c o k i n e t i c s (Marcel Dekker, N e w York, 1982).
275 [2] B.E. Rodda, Analysis of sets of estimates from pharmacokinetic studies, in: Kinetic Data Analysis, ed. L. Endrenyi, pp. 285-298 (Plenum Press, New York, 1981). [3] "C.C. Peck and B.B. Barrett, Nonlinear least-squares regression programs for microcomputers, J. Pharmacokin. Biopharm. 7 (1979) 537-541. [4] Y. Bard, Nonlinear Parameter Estimation (Academic Press, New York, 1974). [5] M. Pfeffer, COMPT, a time sharing program for nonlinear regression analysis of compartmental models of drug distribution, J. Pharmacokin. Biopharm. 1 (1973) 137-163. [6] C.M. Metzler, G.L. Elfring and A.J. McEwen, A package of computer programs for pharmacokinetic modeling, Biometrics 30 (1974) 562. [7] J.G. Wagner, Fundamentals of Clinical Pharmacokinetics (Drug Intelligence, Illinois, 1975). [8] P.V. Pedersen, Curve fitting and modelling in pharmacokinetics and some practical experiences with NONLIN and a new program FUNFIT, J. Pharmacokin. Biopharm, 5 (1977) 513-531. [9] R. Gomeni and C. Gomeni, AUTOMOD: A polyalgorithm for an integrated analysis of linear pharmacokinetic models, Comput. Biol. Med. 9 (1979) 39-48. [10] A. Messori, G.D. Cori and E. Tendi, Iterative least-squares fitting programs in pharmacokinetics for a programmable hand held calculator, Am. J. Hosp. Pharm. 40 (1983) 1673-1684. [11] D.G. Gardner, J.C. Gardner and G. Laush, Method for the analysis of multicomponent exponential decay curves, J. Chem. Phys. 31 (1959) 978-986. [12] M.R. Smith and S.T. Nichols, Improved resolution in the analysis of multicomponent exponential signals, Nuclear lnstrum. Meth. 205 (1983) 479-483.
[13] R.G. Cornell, A method for fitting linear combinations of exponentials, Biometrics 18 (1962) 104-113. [14] D.H. Parsons, Biological problems involving sums of exponential functions of time: an improved method of calculation, Math. Biosci. 9 (1970) 37-47. [15] S.D. Foss, A method of exponential curve fitting by numerical integration, Biometrics 26 (1970) 815-821. [16] K.C. Yeh and K.C. Kwan, Kinetic parameter estimation by numerical algorithms and multiple linear regression: application to pharmacokinetics, J. Pharm. Sci. 68 (1979) 1120-1124. [17] J.R. Koup, Direct linear plotting method for estimation of pharmacokinetics parameters, J. Pharm. Sci. 70 (1981) 1093-1094. [18] A.J. Sedman and J.G. Wagner, CSTRIP, a FORTRAN IV computer program for obtaining polyexponential parameter estimates. J. Pharm. Sci. 65 (1976) 1006-1010. [19] R.D. Brown and J.E. Manno, ESTRIP: a BASIC computer program for obtaining initial polyexponential parameter estimates, J. Pharm. Sci. 67 (1978) 1687-1691. [20] J.G. Leferink and R.A.A. Maes, STRIPACT, an interactive curve fit programme for pharmacokinetic analyses. Arzneim. Forsch. 29 (1979) 1894-1898. [21] P.P. leBlanc and J. Dumas, Calcul des valeurs initiales des parametres pharmacocinetiques a raide d'un calculateur programmable, Therapie 38 (1983) 21-26. [22] A. Dunne and L. Wilson, Multiexponential parameter estimation: elimination of bias, Br. J. Pharmacol. 80 (1983) 714P.