80
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 37. NO. I . JANUARY I Y Y O
A Probabilistic Approach to the Single-Point, Single-Dose Problem WILLIAM L. BRIGGS, ROBERT W . PHELPS,
Absfracf-A general probabilistic approach is applied to the singlepoint, single-dose method for estimating individual infusion rates and serum drug concentrations. By using transformations of probability density functions, the effects of variations in the elimination rate constant upon pharmacokinetic variables may be studied and optimal sampling times may be chosen. Although this study treats the case of error-free sampling in a single-compartment model with a normal distribution of rate constants, the methods presented can be extended to more general situations.
INTRODUCTION HE problem of predicting individual dosage requirements based upon a small number of sampling measurements has received considerable attention in recent years. The use of average or nominal values of pharmacokinetic parameters in predicting dosage rates or steadystate drug concentrations has obvious limitations due to intrapatient variability. In order to incorporate individual variability into such predictions, the “single-point, single-dose’’ was proposed by Slattery et al. [ l ] . In this method, one measurement of serum drug concentration is made at a sampling time after a test dose is given. This single measurement is then used in a linear, one compartment model to predict relevant pharmacokinetic variables. The question that naturally arises is whether an optimal sampling time can be chosen that produces the best possible predictions. Several solutions to this optimal sampling time problem have been given for various predicted variables. Koup [2] discusses the “single-point, single dose” method as it applies to the prediction of steady-state serum drug concentrations. Unadkat and Rowland [3] consider the method for predicting maintenance dosing rates. In a study which most nearly anticipates the approach of our paper, Dossing et al. [4] introduce some statistical structure into the “single-point, single-dose’’ method to determine optimal measurements of clearance. All of these studies determine an optimal sampling time of ts = k,‘ where k, is the average population value of the elimination rate constant. This optimal sampling time is ultimately determined through a linearization argument
T
Manuscript received August 29, 1988; revised April 5 , 1989. W . L. Briggs is with the Department of Mathematics, University of Colorado, Denver, CO 80202. R . W. Phelps and G. D. Swanson are with the Department of Anestheliology, University of Colorado, Denver, CO 80202. IEEE Log Number 8931540.
AND
GEORGE D. SWANSON, MEMBER, IEEE
in which the variation in the predicted variable due to variations in the rate constant is minimized by using the first order terms in a Taylor series expansion. This argument is strictly valid only when variations from k, are small. The work presented in this paper extends the results of these earlier studies by assigning a specific statistical distribution to the elimination rate constant and studying the resulting distribution of various predicted variables. By analyzing these resulting distributions, it is possible to identify optimal sampling times. In certain limiting cases, the optimal sampling time results of previous studies may be recovered. The method does not rely upon a linearization argument or an assumption of small variation. This approach seems to be ideally suited to a variety of pharmacokinetic problems and has applications in other fields as well. Although this study considers only the case of error-free sampling in a single compartment model with a normal distribution of rate constants, the methods presented are being extended to multicompartment models, sampling with error and general distributions of rate constant. PROBLEM FORMULATION A N D SOLUTION We shall consider the situation in which a target drug concentration CT is specified and is to be maintained by a constant infusion rate. (The same question could be asked and the same technique could be applied to the problem of periodic injections.) Using a linear, one-compartment (mono-exponential) model, the serum drug concentration at time t is given by D ~ ( t =) - ( 1 - e - k f ) ,t
vk
>o
where D is the dose rate (mass/time), V is the volume of distribution, and k is the elimination rate constant. Clearly, the steady-state concentration given by this model is CT = D / V k . An initial estimate for the required infusion rate to maintain the target level is simply D = CTVk where, in the absence of any additional information, average population values of V and k would be used. In order to obtain additional, patient-specific information, we allow one measurement of serum drug concentration at time t, following an injection of a test dose d .
0018-9294/90/0100-0080$01.OO
O 1990 IEEE
BRIGGS
PI U / . :
81
APPROACH TO SINGLE-DOSE PROBLEM
In practice, the value of the elimination constant for the ith patient is not known. Therefore, the average value would be used yielding an estimate of the infusion rate as
This measured concentration ( C , ) is given by
This measurement is assumed to be made without error and may be used to determine the volume of distribution V in terms of the other parameters. (The case of measurement with error is important and is under consideration.) An improved estimate of the required infusion rate then follows by solving (2) for V and substituting into the equation for D : (3) At this point, t, is arbitrary and the only value of k available is an average, population value. We now recognize that the elimination constant k has some variability across all patients. This variability in k produces uncertainty in the infusion rate D . Since the sampling time t,vis arbitrary, we ask if there is a choice of t, that minimizes the variation in D due to the variation in k . The variation is minimized when
or when the sampling time t, is given by t, = 1 / k = 1 / k , where k, is the average population value. This argument provides a clear and plausible solution to this optimal sampling problem. However, while the argument accounts for the variation in the rate constant, it does not incorporate a specific distribution of the rate constant across all patients and it does assume the deviation of k from the average population value k, is small. We now assume that the rate constant k has a known distribution and investigate the resulting distribution of the steady-state drug concentrations actually achieved with an optimal dosing strategy. We shall formulate this problem in the following way. Each patient has an individual elimination rate constant and volume of distribution that we denote kj and V i , respectively. If these two quantities were known, then an individual infusion rate Dj = CTk j Vi could be prescribed to reach the target concentration C T . As before, one measurement is allowed following a test dose in order to obtain additional information. Specifically, this measurement of concentration C, taken at sampling time tS can be used to estimate the individual volume Vi as
The resulting serum blood concentration for the ith patient is
where C,yis the measured concentration for the ith patient given by (7)
Combining (6) and (7) gives the concentration ratio
This ratio has the useful property that if a patient is “average” with ki = k,, then r = 1. Furthermore, a deviation away from r = 1 implies an error in the delivered concentration. Minimizing the variance in r minimizes the error in delivered concentrations. At this point, one could impose any sort of distribution upon k . We will assume that k is normally distributed about a known mean k, with a known variance a2; however, other distributions, such as the log-normal, have also been considered. Therefore, the probability density function for k may be written as (9) The problem of finding the resulting density function of r = f ( k ) is standard ( 5 ) . Unfortunately, a solution cannot be given in closed form unless the inverse o f f can be written explicitly. If this cannot be done, as in the present case, then the density function of r must be determined numerically. An algorithm for doing this is given in the Appendix. Even without an explicit representation for the density function of r , it is possible to make some progress toward determining the sampling time t, which minimizes the variance of r . Letting E ( r ) denote the expected value of the random variable r , the variance of r is given by var(r)
~ ( r - ~~ )~ ( r ) .
Furthermore, E ( r 2 )and E ( r ) may be computed knowing only the density function of k since E(r)
This leads to an infusion rate given by
=
=
E(f(k))=
j
m
f(k)4(k)dk --OD
and (4)
E ( r 2 ) = E ( f 2 ( k ) )=
j
m
f2(k)4(k)dk. -m
~
82
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 37, NO. I . JANUARY 1990
Using the definition offgiven in (8) and the assumed density function I$ given in (9), these integrals may be evaluated and simplified. The result is
10.0 - -
8.0- -
1
var ( r ) = 3 e"":{ eU2':(U2+ k: - 4tsk,u2 (10) k, 4 t ; a 4 ) - k: 2t,k,a2 - t , 2 a 4 } .
+
6.0
+
n
.-
4.0 - -
A direct attempt to minimize this expression with respect to ts appears to be futile. However, there are some alternatives that provide valuable insights. In the limit in which o 2 is small, the distribution of rate constants is narrow and the corresponding variance in r may be approximated by var ( r )
U2
7( 1 - 2k,ts k0
r
+ kit:).
Minimization of this expression is straightforward and gives the optimal sampling time of ts = 1 / k , . This agrees with the sampling time obtained earlier without using an assumption about the distribution of rate constants, which also assumes that I k - ko 1 is small. Thus, when the rate constant is distributed narrowly and normally about a mean k,, the sampling time t, = k,' produces the least variance in the concentration ratio. Further information about the distribution of r may be obtained by proceeding numerically. Using the algorithm given in the Appendix, the density function for r = f( k ) , which we denote $( r ) , may be approximated. The result of this calculation is shown in Fig. 1 for the case in which k, = 1 , a 2 = 0 . 0 4 and 0 < ts I2 . As can be seen, the density functions narrow to a sharp peak as ts approaches kc;'. On either side of t,? = k , ' , the distribution of r broadens. All curves along lines of constant t, have unit area. Another perspective on the distribution of r is given in Fig. 2. The full expression for var ( r ) [see ( l o ) ] is evaluated for t, > 0 with U and k, fixed. The value of t, at which var( r ) has a minimum is noted and called t,*. Fig. 2 shows plot of t,* as a function of U for k, = 0.5, k, = 1.0, and k, = 1.5. As can be seen, the optimal sampling time t,* is near k,' when U is small. However, as u 2 increases, the optimal sampling time decreases. We need only remember that the measured concentration taken at time t, is used to determine the relationship between the volume of distribution of a particular patient and the rate constant ki [see ( 2 ) ] . If the measurement is made too early, the value of Vi can be determined more and more exactly, but we know less and less about the value of kj . If the measurement is made too late, the variation of Vi,given various values of k i , increases. Thus, information about ki is lost when ts is too small and information about Vi is lost as the concentration approaches zero. Clearly, there is an optimal sampling time and Fig. 2 suggests that the broader the distribution of rate constants, the sooner the measurement should be made. This
2.
J ts Fig. 1 . Density functions for r = k / k c , e - ' r ( k - k ' where ') k is distributed normally with mean k, = 1 and variance u 2 = 0.04.
20
I .o ko=l.0 k,= 1.5 k,=0.5
I
3
, .
.
.
.
.
.
.
a
'V' -02 .06 . I O .I4 .I8 Fig. 2. The optimal sampling time r,* which minimizes var ( r ) of (9) as a function of U for k , = 0.5, 1.O, 1.5.
is consistent with the concept that with large rate constants, sampling must also be done sooner while concentration levels are still high.
DISCUSSION Previous methods for solving the "single-point, singledose" problem use linearization methods and determine the optimal sampling time by minimizing the sensitivity to changes in the population rate constant. In this paper, we have expanded this approach by assuming a specified distribution for the population rate constants, by determining the resulting distribution for the concentration ratio and then by selecting the optimal sampling time that minimizes the variance of this distribution. The results indicate that when the distribution of k is narrow, there is a high probability that r = 1 when ts = l / k , (Fig. 1).
83
BRIGGS er al.: APPROACH TO SINGLE-DOSE PROBLEM
+
This implies that the delivered concentration Ciis close to the target concentration. This target concentration would be achieved exactly only for the average patient ( ki = ko). These methods appear to be generalizable to a variety of pharmacokinetic problems, many of which are under investigation. The analysis presented in this paper could be applied in a similar way to the problem of determining best estimates of dosing rates for periodic injections. A similar approach is being taken to investigate the effects of sampling error in multiple measurement procedures. The use of a log-normal distribution for the rate constant is also under consideration as well as the use of multicompartment models. This study suggests that the technique is of value whenever one wishes to investigate the distribution of output variables due to a known variation or distribution of noise in input variables.
and grid points y, = g ( a ) n A y for 0 In IN . The number of grid points N determines the accuracy of the approximations to $ and Y. In the calculations of this paper, values of N = 30 and N = 40 were used. By (1 I), we have that
APPENDIX NUMERICAL TRANSFORMATION OF PROBABILITY DENSITY FUNCTIONS Given that a random variable x is distributed according to a density function 4 and a cumulative distribution function X , the problem of finding the resulting density function $ and cumulative distribution function Y of the variable y = g (x) is solved formally in many statistics texts [5]. However, this formal solution must generally be carried out numerically. Let { be differentiable and monotonic increasing on an interval of interest a Ix Ib. Also let P ( y ' Iy ) denote the probability that y ' Iy . By the definition of the cumulative distribution function, we have that
Therefore, combining (14) and (15) leads to
Y( Y )
P ( Y'
IY
=
'("
I -I(
=
"(')
X"
Xn- I
=
Y ( Y f l - I )+
j
444 dt
Xn
4 ( t ) dt. Xn- I
If the interval a Ix Ib is chosen such that 4 ( a ) is small, then it is reasonable to take Y ( y o ) = Y [ ( a ) ] = 0. With Y ( y o ) = 0 , (16) may be used to successively compute Y ( y l ) , * * . , Y ( y N ) , as long as the integral in (16) can be approximated. The simplest way to approximate the integral is by the trapezoid rule:
) = P ( g (x') 5 Y )
=
y ) ) = x( g
-I(
y)).
')
Since the cumulative distribution function and the probability density function are related by
'(')
Equation (12) applied to 4 and X gives
Or
'(')
=
s'
-m
'(t)dt'
(12)
we have
The density function 4 is known. However, the coordinates xn are not known and must be computed. These COordinates are the roots of the equation g (x,) = Y,, where the grid points yn are known. This root finding problem is most easily solved by Newton's method with an initial guess o f x n - l . Therefore, the sequence Y ( x , ) f o r n = 1, 2, * ' . , N may be computed by (16) in the form
YfYn) = U Y n - I ) + +(4bn,
(17)
+ 4(xn-d)(xn - 4 - 1 ) . =
X
(g
-I
(y))
d 6 (6
( y )).
( 13)
If the function g can be inverted explicitly, then it may be possible to compute $ and Y as given by (1 1) and (13). More likely, 6 cannot be inverted analytically and the following numerical scheme is one of many possible methods to approximate $ and Y . A uniform grid is placed on the interval { ( a ) Iy I 1( b )with grid spacing
Ay =
- {(a)
N
This determines the cumulative distribution function for y = g (x). Using (12) and (17), we may approximate the density function $ by $(Yn-I/d
=
Y(yri) - Y ( y n - 1 1 AY (xn - xrz-l)
(18)
for I In I N . Equations (17) and (18) generate approximations to Y and $ on the chosen grid { yo. y l , . . . ,
84
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 37, NO. I. JANUARY 1990
y N }. The procedure may be summarized in the following algorithm:
Initialize: A y = g ( b ) N
g(4
g ( a ) + nAy,
Y ( Y o ) = 0, xo = a , F o r n = 1,2, Solve g (x,) guess X n - I
-
e
=
mathematics at the University of Colorado in 1971 and the Ph.D. degree in applied mathematics at Harvard University, Cambridge, MA, in 1976. His research interests include parallel computation and mathematical problems in the biosciences. He is currently in the Department of Mathematics, University of Colorado at Denver.
g, 4.
Specify: a , b, N ,
yn =
William L. Briggs received the B.A. degree in
*
0 d+O)
In IN
=0
Robert W. Phelps received the Ph.D. degree in biophysics from Stanford University, Stanford, CA, in 1972, and the M.D. degree from the University of Rochester, Rochester, NY, in 1979. He joined the faculty of the University of Colorado Health Sciences Center in 1982 and is now an Associate Professor of Anesthesiology. Although he is currently engaged mostly in teaching all types of clinical anesthesia, his interests include anesthesia monitoring equipment, computer applications and mathematical simulation of clin-
,N yn by Newton’s method with initial
-1)
- x n -- 1 ) .
If g is monotone decreasing on a Ix Ib, the algorithm may be used with minor modifications. If g is not monotone on the interval of interest (as is the case in the problem of this paper), the interval must be divided into subintervals on which 6 is monotone. The results of these subproblems may then be combined to produce Y and 4.
REFERENCES
ical phenomena
George D. Swanson (S’69-M’72) was born in Eureka, CA, on September 19, 1942. He received the B.S. degree from California State Polytechnic University, San Luis Obispo, in 1965, and the M.S. and Ph.D. degrees from Stanford, CA in 1966 and 1972, respectively. From 1968 to 1.969 he was an Acting Instructor in the Department of Electrical Engineering, Stanford University, and from 1968 to 1972 he was a Research AssistantiAssociate in the DeDartment of Anesthesioloev. Stanford Universitv. In 1972, he joined the Dkpartment of Anesthesiology and the Department of System Science, University of California, Los Angeles, as an Assistant Professor, where he helped develop a research and teaching program in the application of system science to medicine and physiology. In 1975 he was Visiting Assistant Professor of Anesthesiology, University of Colorado Medical School, Denver. In 1976 he joined the faculty of the Department of Anesthesiology and the Department of Biometrics, University of Colorado Health Sciences Center, Denver, where he is an Associate Professor and Director of Anesthesiology Research. Most recently, he has been appointed professor of Exercise Physiology at California State University, Chico. Dr. Swanson is a member of the American Physiological Society, the Sigma Xi Research Society, the Bioengineering Society, the American College of Sports Medicine, and the Aerospace Medical Association. -I
[I] J . T. Slattery, M. Gibaldi, and J . R. Koup, “Prediction of drug concentration at steady state from a single determination of concentration after an initial dose,” Clin. Pharmacokinefics, vol. 5 , pp. 377-385, 1980. [2] J . R. Koup, “Single-point prediction methods: A critical review,” Drug Intell. Clin. Pharma., vol. 16, pp. 855-862, 1982. [3] J . D. Unadkat and M. Rowland, “Further considerations of the singlepoint single-dose method to estimate individual maintenance dosage requirements,” Ther. Drug Monit, vol. 4, pp. 201-208, 1982. [4] M. Dossing, A. Volund, and H. E. Poulsen, “Optimal sampling times for minimum variance of clearance determination,” Brit. J . Clin. Pharrnacol., vol. 15, pp. 231-235, 1983. [5] R. V . Hogg and A. T. Craig, Introduction to Mathematical Statisrics, 4th ed. New York: MacMillan, 1978.
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING. VOL. 37. NO. I . JANUARY I Y Y O
A Probabilistic Approach to the Single-Point, Single-Dose Problem WILLIAM L. BRIGGS, ROBERT W . PHELPS,
Absfracf-A general probabilistic approach is applied to the singlepoint, single-dose method for estimating individual infusion rates and serum drug concentrations. By using transformations of probability density functions, the effects of variations in the elimination rate constant upon pharmacokinetic variables may be studied and optimal sampling times may be chosen. Although this study treats the case of error-free sampling in a single-compartment model with a normal distribution of rate constants, the methods presented can be extended to more general situations.
INTRODUCTION HE problem of predicting individual dosage requirements based upon a small number of sampling measurements has received considerable attention in recent years. The use of average or nominal values of pharmacokinetic parameters in predicting dosage rates or steadystate drug concentrations has obvious limitations due to intrapatient variability. In order to incorporate individual variability into such predictions, the “single-point, single-dose’’ was proposed by Slattery et al. [ l ] . In this method, one measurement of serum drug concentration is made at a sampling time after a test dose is given. This single measurement is then used in a linear, one compartment model to predict relevant pharmacokinetic variables. The question that naturally arises is whether an optimal sampling time can be chosen that produces the best possible predictions. Several solutions to this optimal sampling time problem have been given for various predicted variables. Koup [2] discusses the “single-point, single dose” method as it applies to the prediction of steady-state serum drug concentrations. Unadkat and Rowland [3] consider the method for predicting maintenance dosing rates. In a study which most nearly anticipates the approach of our paper, Dossing et al. [4] introduce some statistical structure into the “single-point, single-dose’’ method to determine optimal measurements of clearance. All of these studies determine an optimal sampling time of ts = k,‘ where k, is the average population value of the elimination rate constant. This optimal sampling time is ultimately determined through a linearization argument
T
Manuscript received August 29, 1988; revised April 5 , 1989. W . L. Briggs is with the Department of Mathematics, University of Colorado, Denver, CO 80202. R . W. Phelps and G. D. Swanson are with the Department of Anestheliology, University of Colorado, Denver, CO 80202. IEEE Log Number 8931540.
AND
GEORGE D. SWANSON, MEMBER, IEEE
in which the variation in the predicted variable due to variations in the rate constant is minimized by using the first order terms in a Taylor series expansion. This argument is strictly valid only when variations from k, are small. The work presented in this paper extends the results of these earlier studies by assigning a specific statistical distribution to the elimination rate constant and studying the resulting distribution of various predicted variables. By analyzing these resulting distributions, it is possible to identify optimal sampling times. In certain limiting cases, the optimal sampling time results of previous studies may be recovered. The method does not rely upon a linearization argument or an assumption of small variation. This approach seems to be ideally suited to a variety of pharmacokinetic problems and has applications in other fields as well. Although this study considers only the case of error-free sampling in a single compartment model with a normal distribution of rate constants, the methods presented are being extended to multicompartment models, sampling with error and general distributions of rate constant. PROBLEM FORMULATION A N D SOLUTION We shall consider the situation in which a target drug concentration CT is specified and is to be maintained by a constant infusion rate. (The same question could be asked and the same technique could be applied to the problem of periodic injections.) Using a linear, one-compartment (mono-exponential) model, the serum drug concentration at time t is given by D ~ ( t =) - ( 1 - e - k f ) ,t
vk
>o
where D is the dose rate (mass/time), V is the volume of distribution, and k is the elimination rate constant. Clearly, the steady-state concentration given by this model is CT = D / V k . An initial estimate for the required infusion rate to maintain the target level is simply D = CTVk where, in the absence of any additional information, average population values of V and k would be used. In order to obtain additional, patient-specific information, we allow one measurement of serum drug concentration at time t, following an injection of a test dose d .
0018-9294/90/0100-0080$01.OO
O 1990 IEEE
BRIGGS
PI U / . :
81
APPROACH TO SINGLE-DOSE PROBLEM
In practice, the value of the elimination constant for the ith patient is not known. Therefore, the average value would be used yielding an estimate of the infusion rate as
This measured concentration ( C , ) is given by
This measurement is assumed to be made without error and may be used to determine the volume of distribution V in terms of the other parameters. (The case of measurement with error is important and is under consideration.) An improved estimate of the required infusion rate then follows by solving (2) for V and substituting into the equation for D : (3) At this point, t, is arbitrary and the only value of k available is an average, population value. We now recognize that the elimination constant k has some variability across all patients. This variability in k produces uncertainty in the infusion rate D . Since the sampling time t,vis arbitrary, we ask if there is a choice of t, that minimizes the variation in D due to the variation in k . The variation is minimized when
or when the sampling time t, is given by t, = 1 / k = 1 / k , where k, is the average population value. This argument provides a clear and plausible solution to this optimal sampling problem. However, while the argument accounts for the variation in the rate constant, it does not incorporate a specific distribution of the rate constant across all patients and it does assume the deviation of k from the average population value k, is small. We now assume that the rate constant k has a known distribution and investigate the resulting distribution of the steady-state drug concentrations actually achieved with an optimal dosing strategy. We shall formulate this problem in the following way. Each patient has an individual elimination rate constant and volume of distribution that we denote kj and V i , respectively. If these two quantities were known, then an individual infusion rate Dj = CTk j Vi could be prescribed to reach the target concentration C T . As before, one measurement is allowed following a test dose in order to obtain additional information. Specifically, this measurement of concentration C, taken at sampling time tS can be used to estimate the individual volume Vi as
The resulting serum blood concentration for the ith patient is
where C,yis the measured concentration for the ith patient given by (7)
Combining (6) and (7) gives the concentration ratio
This ratio has the useful property that if a patient is “average” with ki = k,, then r = 1. Furthermore, a deviation away from r = 1 implies an error in the delivered concentration. Minimizing the variance in r minimizes the error in delivered concentrations. At this point, one could impose any sort of distribution upon k . We will assume that k is normally distributed about a known mean k, with a known variance a2; however, other distributions, such as the log-normal, have also been considered. Therefore, the probability density function for k may be written as (9) The problem of finding the resulting density function of r = f ( k ) is standard ( 5 ) . Unfortunately, a solution cannot be given in closed form unless the inverse o f f can be written explicitly. If this cannot be done, as in the present case, then the density function of r must be determined numerically. An algorithm for doing this is given in the Appendix. Even without an explicit representation for the density function of r , it is possible to make some progress toward determining the sampling time t, which minimizes the variance of r . Letting E ( r ) denote the expected value of the random variable r , the variance of r is given by var(r)
~ ( r - ~~ )~ ( r ) .
Furthermore, E ( r 2 )and E ( r ) may be computed knowing only the density function of k since E(r)
This leads to an infusion rate given by
=
=
E(f(k))=
j
m
f(k)4(k)dk --OD
and (4)
E ( r 2 ) = E ( f 2 ( k ) )=
j
m
f2(k)4(k)dk. -m
~
82
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 37, NO. I . JANUARY 1990
Using the definition offgiven in (8) and the assumed density function I$ given in (9), these integrals may be evaluated and simplified. The result is
10.0 - -
8.0- -
1
var ( r ) = 3 e"":{ eU2':(U2+ k: - 4tsk,u2 (10) k, 4 t ; a 4 ) - k: 2t,k,a2 - t , 2 a 4 } .
+
6.0
+
n
.-
4.0 - -
A direct attempt to minimize this expression with respect to ts appears to be futile. However, there are some alternatives that provide valuable insights. In the limit in which o 2 is small, the distribution of rate constants is narrow and the corresponding variance in r may be approximated by var ( r )
U2
7( 1 - 2k,ts k0
r
+ kit:).
Minimization of this expression is straightforward and gives the optimal sampling time of ts = 1 / k , . This agrees with the sampling time obtained earlier without using an assumption about the distribution of rate constants, which also assumes that I k - ko 1 is small. Thus, when the rate constant is distributed narrowly and normally about a mean k,, the sampling time t, = k,' produces the least variance in the concentration ratio. Further information about the distribution of r may be obtained by proceeding numerically. Using the algorithm given in the Appendix, the density function for r = f( k ) , which we denote $( r ) , may be approximated. The result of this calculation is shown in Fig. 1 for the case in which k, = 1 , a 2 = 0 . 0 4 and 0 < ts I2 . As can be seen, the density functions narrow to a sharp peak as ts approaches kc;'. On either side of t,? = k , ' , the distribution of r broadens. All curves along lines of constant t, have unit area. Another perspective on the distribution of r is given in Fig. 2. The full expression for var ( r ) [see ( l o ) ] is evaluated for t, > 0 with U and k, fixed. The value of t, at which var( r ) has a minimum is noted and called t,*. Fig. 2 shows plot of t,* as a function of U for k, = 0.5, k, = 1.0, and k, = 1.5. As can be seen, the optimal sampling time t,* is near k,' when U is small. However, as u 2 increases, the optimal sampling time decreases. We need only remember that the measured concentration taken at time t, is used to determine the relationship between the volume of distribution of a particular patient and the rate constant ki [see ( 2 ) ] . If the measurement is made too early, the value of Vi can be determined more and more exactly, but we know less and less about the value of kj . If the measurement is made too late, the variation of Vi,given various values of k i , increases. Thus, information about ki is lost when ts is too small and information about Vi is lost as the concentration approaches zero. Clearly, there is an optimal sampling time and Fig. 2 suggests that the broader the distribution of rate constants, the sooner the measurement should be made. This
2.
J ts Fig. 1 . Density functions for r = k / k c , e - ' r ( k - k ' where ') k is distributed normally with mean k, = 1 and variance u 2 = 0.04.
20
I .o ko=l.0 k,= 1.5 k,=0.5
I
3
, .
.
.
.
.
.
.
a
'V' -02 .06 . I O .I4 .I8 Fig. 2. The optimal sampling time r,* which minimizes var ( r ) of (9) as a function of U for k , = 0.5, 1.O, 1.5.
is consistent with the concept that with large rate constants, sampling must also be done sooner while concentration levels are still high.
DISCUSSION Previous methods for solving the "single-point, singledose" problem use linearization methods and determine the optimal sampling time by minimizing the sensitivity to changes in the population rate constant. In this paper, we have expanded this approach by assuming a specified distribution for the population rate constants, by determining the resulting distribution for the concentration ratio and then by selecting the optimal sampling time that minimizes the variance of this distribution. The results indicate that when the distribution of k is narrow, there is a high probability that r = 1 when ts = l / k , (Fig. 1).
83
BRIGGS er al.: APPROACH TO SINGLE-DOSE PROBLEM
+
This implies that the delivered concentration Ciis close to the target concentration. This target concentration would be achieved exactly only for the average patient ( ki = ko). These methods appear to be generalizable to a variety of pharmacokinetic problems, many of which are under investigation. The analysis presented in this paper could be applied in a similar way to the problem of determining best estimates of dosing rates for periodic injections. A similar approach is being taken to investigate the effects of sampling error in multiple measurement procedures. The use of a log-normal distribution for the rate constant is also under consideration as well as the use of multicompartment models. This study suggests that the technique is of value whenever one wishes to investigate the distribution of output variables due to a known variation or distribution of noise in input variables.
and grid points y, = g ( a ) n A y for 0 In IN . The number of grid points N determines the accuracy of the approximations to $ and Y. In the calculations of this paper, values of N = 30 and N = 40 were used. By (1 I), we have that
APPENDIX NUMERICAL TRANSFORMATION OF PROBABILITY DENSITY FUNCTIONS Given that a random variable x is distributed according to a density function 4 and a cumulative distribution function X , the problem of finding the resulting density function $ and cumulative distribution function Y of the variable y = g (x) is solved formally in many statistics texts [5]. However, this formal solution must generally be carried out numerically. Let { be differentiable and monotonic increasing on an interval of interest a Ix Ib. Also let P ( y ' Iy ) denote the probability that y ' Iy . By the definition of the cumulative distribution function, we have that
Therefore, combining (14) and (15) leads to
Y( Y )
P ( Y'
IY
=
'("
I -I(
=
"(')
X"
Xn- I
=
Y ( Y f l - I )+
j
444 dt
Xn
4 ( t ) dt. Xn- I
If the interval a Ix Ib is chosen such that 4 ( a ) is small, then it is reasonable to take Y ( y o ) = Y [ ( a ) ] = 0. With Y ( y o ) = 0 , (16) may be used to successively compute Y ( y l ) , * * . , Y ( y N ) , as long as the integral in (16) can be approximated. The simplest way to approximate the integral is by the trapezoid rule:
) = P ( g (x') 5 Y )
=
y ) ) = x( g
-I(
y)).
')
Since the cumulative distribution function and the probability density function are related by
'(')
Equation (12) applied to 4 and X gives
Or
'(')
=
s'
-m
'(t)dt'
(12)
we have
The density function 4 is known. However, the coordinates xn are not known and must be computed. These COordinates are the roots of the equation g (x,) = Y,, where the grid points yn are known. This root finding problem is most easily solved by Newton's method with an initial guess o f x n - l . Therefore, the sequence Y ( x , ) f o r n = 1, 2, * ' . , N may be computed by (16) in the form
YfYn) = U Y n - I ) + +(4bn,
(17)
+ 4(xn-d)(xn - 4 - 1 ) . =
X
(g
-I
(y))
d 6 (6
( y )).
( 13)
If the function g can be inverted explicitly, then it may be possible to compute $ and Y as given by (1 1) and (13). More likely, 6 cannot be inverted analytically and the following numerical scheme is one of many possible methods to approximate $ and Y . A uniform grid is placed on the interval { ( a ) Iy I 1( b )with grid spacing
Ay =
- {(a)
N
This determines the cumulative distribution function for y = g (x). Using (12) and (17), we may approximate the density function $ by $(Yn-I/d
=
Y(yri) - Y ( y n - 1 1 AY (xn - xrz-l)
(18)
for I In I N . Equations (17) and (18) generate approximations to Y and $ on the chosen grid { yo. y l , . . . ,
84
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 37, NO. I. JANUARY 1990
y N }. The procedure may be summarized in the following algorithm:
Initialize: A y = g ( b ) N
g(4
g ( a ) + nAy,
Y ( Y o ) = 0, xo = a , F o r n = 1,2, Solve g (x,) guess X n - I
-
e
=
mathematics at the University of Colorado in 1971 and the Ph.D. degree in applied mathematics at Harvard University, Cambridge, MA, in 1976. His research interests include parallel computation and mathematical problems in the biosciences. He is currently in the Department of Mathematics, University of Colorado at Denver.
g, 4.
Specify: a , b, N ,
yn =
William L. Briggs received the B.A. degree in
*
0 d+O)
In IN
=0
Robert W. Phelps received the Ph.D. degree in biophysics from Stanford University, Stanford, CA, in 1972, and the M.D. degree from the University of Rochester, Rochester, NY, in 1979. He joined the faculty of the University of Colorado Health Sciences Center in 1982 and is now an Associate Professor of Anesthesiology. Although he is currently engaged mostly in teaching all types of clinical anesthesia, his interests include anesthesia monitoring equipment, computer applications and mathematical simulation of clin-
,N yn by Newton’s method with initial
-1)
- x n -- 1 ) .
If g is monotone decreasing on a Ix Ib, the algorithm may be used with minor modifications. If g is not monotone on the interval of interest (as is the case in the problem of this paper), the interval must be divided into subintervals on which 6 is monotone. The results of these subproblems may then be combined to produce Y and 4.
REFERENCES
ical phenomena
George D. Swanson (S’69-M’72) was born in Eureka, CA, on September 19, 1942. He received the B.S. degree from California State Polytechnic University, San Luis Obispo, in 1965, and the M.S. and Ph.D. degrees from Stanford, CA in 1966 and 1972, respectively. From 1968 to 1.969 he was an Acting Instructor in the Department of Electrical Engineering, Stanford University, and from 1968 to 1972 he was a Research AssistantiAssociate in the DeDartment of Anesthesioloev. Stanford Universitv. In 1972, he joined the Dkpartment of Anesthesiology and the Department of System Science, University of California, Los Angeles, as an Assistant Professor, where he helped develop a research and teaching program in the application of system science to medicine and physiology. In 1975 he was Visiting Assistant Professor of Anesthesiology, University of Colorado Medical School, Denver. In 1976 he joined the faculty of the Department of Anesthesiology and the Department of Biometrics, University of Colorado Health Sciences Center, Denver, where he is an Associate Professor and Director of Anesthesiology Research. Most recently, he has been appointed professor of Exercise Physiology at California State University, Chico. Dr. Swanson is a member of the American Physiological Society, the Sigma Xi Research Society, the Bioengineering Society, the American College of Sports Medicine, and the Aerospace Medical Association. -I
[I] J . T. Slattery, M. Gibaldi, and J . R. Koup, “Prediction of drug concentration at steady state from a single determination of concentration after an initial dose,” Clin. Pharmacokinefics, vol. 5 , pp. 377-385, 1980. [2] J . R. Koup, “Single-point prediction methods: A critical review,” Drug Intell. Clin. Pharma., vol. 16, pp. 855-862, 1982. [3] J . D. Unadkat and M. Rowland, “Further considerations of the singlepoint single-dose method to estimate individual maintenance dosage requirements,” Ther. Drug Monit, vol. 4, pp. 201-208, 1982. [4] M. Dossing, A. Volund, and H. E. Poulsen, “Optimal sampling times for minimum variance of clearance determination,” Brit. J . Clin. Pharrnacol., vol. 15, pp. 231-235, 1983. [5] R. V . Hogg and A. T. Craig, Introduction to Mathematical Statisrics, 4th ed. New York: MacMillan, 1978.
