- Email: [email protected]
Segmentally Continuous input Functions in Linear Multicompartment Systems
Segmentally Continuous Input Functions in Linear Multicompartment Systems RICHARD SUVERKRUP Received November 21, 1983, from the Pharmazeufisches Insrirur, Phairrnazeutische Technolcgie der Rheinischen Friedrich-Wilhelms-Universit&, 5300 Bonn-Endenich, F.R.G. Accepted for publication July 6,1984. Abstract 0 The convolutions of segmentally continuous zero- and firstorder input functions with the general form of pharmacokinetic multicom-
partmental disposition functions can be used to analyze and simulate the time course of drug invasion into the central compartment of mammiliary models and certain physiologically relevant recirculating systems. The generalized model equations may be used to assess the reliability and intercorrelations of parameter estimates directly, since partial derivatives with respect to all model constants can be computed explicitly. In combination with curve-fitting algorithms, input functions identical to those of the point-area deconvolution can be obtained, although at the expense of significantly longer computing times. On the other hand, the range of potential applications goes far beyond the reconstruction of the time course of drug absorption. Smooth input functions are, in many cases, only a rough approximation of rather erratic absorption profiles; the time course of absorption is possibly subject to greater intra- and intersubject variability than most other processes of pharmacokinetic interest, but it is a prerequisite for the rationai design of dosage forms. The algorithms commonly employed to analyze drug invasion characteristics can be grouped into three classes: 1. Curve fitting of explicit convolutions of simple input functions with the disposition functions of compartmental models, which may be specified in the form of microscopic rate constants or as parameters of the general solution for plasma concentration profiles after impulse input, represent the first class. Typical examples are the functions given by Teorell,'~2 C ~ t l e r ,and ~ . ~Veng-Peder~en.~-' Generally suitable input functions should be flexible enough to describe the general characteristics of drug invasion in sufficient detail, yet so stiff that single measurements which do not fit into the overall pattern do not cause them to go astray. 2. Algorithms based on mass-balance considerations like the methods of Wagner and Nelson9 and Loo and Riegelman" are model independent with respect to input kinetics but rely on assumptions concerning the structure of the disposition model, since microscopic rate constants are used in the course of computations. Vaughan and Dennis" have demonstrated that the Loo-Riegelman method remains operative even if these assumptions do not strictly hold. 3. Numerical deconvolution procedures like the point-area method of Vaughan and Dennis," midpoint-impulse of Chiou,I3 and the trapezoidal approximation of Stepanek,I4 which is practicable even if a computer is not available, are derived from linear systems theory and do not imply any prerequisites with respect to the type of input kinetics nor specific properties of the model beyond its linearity and time independence. Experimental error is attributed directly to the input function and may give rise to temporarily negative apparent absorption rates in the vicinity of the peak, at the plateau, or on the downslope of the plasma level curve, but not while plasma levels rise markedly, unless error is large. 136
1 Journal of Pharmaceutical Sciences Vot. 74, No. 2, February 1985
The problem is aggravated if plasma levels after intravenous bolus injection are used directly as the weighting function. This is possible if equally spaced samples have been taken at the same times after intravenous administration and absorption. Under these conditions, most deconvolution algorithms assume a particularly simple and attractive form. The shortcut suggests itself in the case of the trapezoidal approximation, which is equivalent to the point-area method if fractional areas under the weighting function are computed by the trapezoidal rule. If this is done, the errors of two data sets are compounded and the outcome is uncertain. Sampling schedules employed in practice rarely coincide after intravenous administration and absorption, and samples are hardly ever taken at equal intervals. Therefore, in order to minimize the computing effort, data are frequently interpolated. At the same time, both the impulse response and the response function to be deconvoluted may be smoothed to make the input profile look more acceptable. This appears to be largely a cosmetic operation, which should be considered carefully, since systematic errors may be introduced even if the data are flawless.
Theoretical Section A general disposition function for linear mammillary models was given by Benet" and extended by Vaughan and cow o r k e r ~ ' ~to . ' ~include input into peripheral compartments. The matrix approach to Laplace transform solutions for linear multicompartment systems arbitrarily interconnected by firstorder transfer has been developed by Wolf, Heinzel, and coworker~'~.'~ and by Veng-Pedersen?' whose notation will be used below. Although, in the general case, input may occur at any time into any compartment with variable rates, only the particular situation where input is restricted to one compartment and enters the system as an impulse at time zero will be considered here, in order to define the weighting function of the system with respect to this compartment. The differential equations describing disposition in a dense n-compartment system with rate constants greater or equal to zero may be written:
+
n
c + ... +
+
j=l j#i
+ ...
+ kin
... -
n-1
c
j=1
0022-3549/85/0200-0136$0 1 .0010 0 1985, American PharmaceuticalAssociation
or briefly: x ' = (KT- Z)x
(2)
K Tis the n x n matrix with rate constants kGfor drug transfer from the ith to the j t h compartment, i = 1 , 2 , . . + ,n and j = 1, 2, . . . , n, where k, = 0 for i = j and Z is a diagonal matrix with elements Ej: n
E.= C k..
compartment model during and after infusion with rate ri, may be written:
In the case of first-order absorption, a lag time to may be taken into account, by which the continuous input response function is also shifted in time, and the first-step variable is defined as:
i=l i#j
With initial conditions xl.o = D and xi.0 = (i = 2, 3, . . .,n),the solution for x1 in the Laplace domain can be written:
D -fl(s)=
(s
+ Ej)
j-2 n
n (s +
In the definition of {I, a scale reduction due to the removal of drug from the absorption site in the time interval to+ has to be introduced, and we have:
(3) Ah)
k-1
where X k are the n eigenvalues of - Z. Dividing the inverse transform by the scaling factor V , , which relates amounts in this compartment to concentrations, the unit-impulse response function is:
Thus, the time course of drug concentrations in plasma for the one-compartment open body model with truncated first-order absorption is:
n
k=l
k#i
In pharmacokinetics, further restrictions with respect to the type of models are common: complex eigenvalues may arise if compartments are connected in loops, resulting in damped sinusoidal oscillations in open systems. Since this is hardly observed in reality, interest has been focused on mammillary models in the broad sense, where branched systems of peripheral compartments may be connected to the central compartment but lateral interconnections forming loops are absent.'l A solution of the form given above can always be derived starting from a model of given structure and known microscopic rate constants, but it is not possible to define both structure and values of rate constants unambiguously from the parameters E,,X i , and V , of the general solution. The problem of vanishing exponentials has been discussed by Wagner22and by Vaughan and who have shown that the number of models corresponding to any given system function is countably infinite if only plasma concentrations are measured. This family of models may, however, be mapped uniquely onto its simplest member, the n-compartment mammillary model with elimination only from the central compartment. Input Function-If a zero- or first-order drug input into a linear compartmental system exhibiting multiexponential input response is stopped abruptly, the time course of plasma concentrations can be represented as the sum of a continuous input response function and a truncation correction function of identical form.24In the case of zero-order infusion starting at time zero, the only differences are the sign of the corrective function and its displacement along the time axis by t l , the time when infusion stops. The notation becomes particularly Let 5; be simple if two-value step variables are intr~duced.'~ defined as: (5)
Then, the time course of plasma concentrations in the one-
If the invasion rate (rin)or the absorption rate constant (k.) do not drop to zero at time t l , but assume values of rin.2or k,.2,respectively, differing from the corresponding values of the first phase (which are now designated and k,.l), a second phase follows. The times at which rate constants are assumed to change their values may be called transition times. For zero-order invasion, a second step variable has to be defined, which changes its value at time t2:
Then, eq. 6 can be amended to read
Similarly, for first-order absorption:
For the case of m-phasic zero-order invasion, a general input function in the domain of the Laplace transform can be derived
The corresponding function for m-phasic first-order absorption is somewhat more complex, but still transparent:
x
[e-~ta-l
-
e-ka,(t.-t.-l)
e
-3t,
(14) l
The convolution product of these segmentally continuous Journal of Pharmaceutical Sciences Vol. 74. No. 2, February 1985
/ 137
input functions with m invasion phases and the general disposition function for n-compartment systems is readily farmed. After inverse transformation into the time domain, the plasma level function for multiphasic zero-order invasion reads:
i=l
- ti
with: (16)
Although the function looks unwieldy a t first glance, its structure is highly regular. The zero-order input response function consists of a time-independent steady-state term and n negative exponential terms. Both numerators and denominators of their coefficients are continued products formed after the same pattern. Note that the steady-state term may be converted to the same form by introducing a virtual constant Xo = 0. For each complete invasion phase, there is one continuous input and one corrective partial function. The corresponding plasma level function for segmentally continuous first-order input reads:
derivatives with respect to all model parameters can be computed explicitly. Enterohepatic Recirculation-The approach developed for mammillary models in the broad sense is also applicable to certain types of models with closed loops describing enterohepatic recirculation. Although the production of bile is a continuous process, bile flow in humans is not steady, as in some animal species. In humans, the emptying of the gallbladder is triggered by stimuli related to food intake. In the context of multicompartmental systems, the contents of the gallbladder may be considered external to the body, like the contents of the urinary bladder or of the gastrointestinal tract. In order to keep the following example simple, first-order transfers are assumed throughout (for the time being), with the only exception of bile release into the duodenum. Then, with one-compartment disposition kinetics and assuming the liver to be a compartment of vanishing volume, Fig. 1is a schematic representation of the model. As long as bile accumulates in the gallbladder, the switch between compartments 4 and 0 is open, and the following system of differential equations holds prior to the first emptying at time tb x; =
-(ko1
+ kz + ko4)Xo
x: =
kolX0
- (kl2 + k13 + k14)XI
x; =
kol xo
hzXi
xi =
ki3Xi
kQ4xo
x; =
(19)
If a dose D is introduced into compartment 1 as a bolus, we have initial conditions a t t = 0, x1 = D,xo = x2 = x3 = x4 o = 0. Immediately prior to the emptying at time tb the distribution of drug amounts in the system is:
-
n
x3
= ki3
x4 = k14
+E
kl4xl
1'"' 1"'
x1 dt
x1 dt
(20)
For the sake of simplicity, again it is assumed that the whole amount of bile accumulated up to t b . 1 is released in so short a time interval ( A t ) that the change of amounts of drugs in compartments 1,2,and 3 is negligible. Subsequently, the transport from compartment 4 to compartment 0 is inhibited and eq. 19 can be solved until the gallbladder is emptied the next e
(17)
1*1
with:
(18)
Because of their repetitive pattern of formation, these equations are readily implemented in the form of compact, yet structurally transparent, computer programs. Even their partial 138
Journal of Pharmaceutical Sciences Vol. 74, No. 2, February 1985
Figure 1-Compartmental model for enterohepatic recirculation. Key to compartments: (0) GI tract; (1) general circulation; (2) metaboiites; (3) drug excreted unchanged into urine; (4) gallbladder.
Table I-Comparison of Staircase Invasion Functions Obtained by Point-Area Deconvolution and Nonlinear Regression Using a Multiphasic Zero-Order Invasion Model with True Values
time at tb.2. Now, the initial conditions are:
0.0 1 .o
2.0 3.0 4.0 5.0 6.0 7.0 X4.tbl
(21)
=0
It is convenient to define the collective rate constants for drug elimination as k, = klz + kI3 + k14 and for its removal from the GI tract as k. = kl+ k2+ k4. The solution for x1 is obvious; superimposed upon the primary impulse response function: = D.e-%'
(22)
and with f = k,,,/k., we have for the amount of recirculated drug after the first release of bile:
0.0 4.25 5.31 5.06 4.37 3.60 2.89 2.30
0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0078
-
0.4912 0.2463 0.1230 0.0622 0.0311 0.0147 0.0077
-
0.4912 0.2463 0.1230 0.0622 0.0306 0.0155 0.0075
-
computer-generated data with and without random error added and from experimental data, taking the times of observations as transition times. The staircase input functions arrived at were virtually identical with those obtained by the point-area deconvolution method of Vaughan and Dennis. The unit-impulse response function given in their paper:"
G ( t ) = 5.Oe-O.*' + 5.0e-'.'
(24)
can be formulated according to eq. 4 as:
Its convolution with the first-order input function: inl(t) = 0.693e-0693' The total amount of drug in compartment 1 is XI = x1 , + x1 r ~ . From time t b on, the intermediate reservoir in compartment 4 is fed from two sources: besides the biliary accumulation of drug cleared from the general circulation, a first-pass process may contribute significantly. Looking only at compartment 1 (drug in plasma) and possibly at compartment 3 (unchanged drug excreted into urine), the system behaves as if a second dose D, has been introduced into compartment 0 at time tb 1. The extension to more than one instance of bile release and the multicompartmetal models of drug disposition is conceptually simple. Several simplifications have been introduced into the example for the sake of clarity, which may be dropped in the general case. Since compartment 0 is considered to be external and a predecessor of compartment 1, the linearity of the model is maintained, as far as distribution and elimination are concerned, if first-order processes of drug absorption, presystemic metabolism, and presystemic biliary secretion are replaced by arbitrary transfer functions, all greater than or equal to zero at any time. Moreover, instantaneous and complete emptying of the gallbladder is not a strict requirement if interest is focused on the time course of absorption of recirculated drug into compartment 1. In order to define the disposition characteristics of a multicompartment system adequately after intravenous administration, i.e., to estimate the parameters E, and h, of its general solution for its circulation compartment with reasonable precision, there must only be sufficiently long periods in the initial and terminal phases during which release of bile and reabsorption are negligible.
Experimental Section For the purpose of data analysis, a subroutine incorporating eqs. 15-18 has been combined with the Nelder-Mead simplex algorithm for minimization of the sum of (weighted) residual squares.26This method, although less efficient than other optimization procedures, has been chosen because it allows the user to fix the values of known model constants conveniently by specifying zero initial variation steps for these parameters. Comparison wit h Point-Area Deconvolution-Segmentally continuous zero-order input rates were estimated from
(26)
was analyzed both by deconvolution and curve fitting, holding the values of V ] , EP,XI, and hz constant. The zero-order invasion rates obtained. for each interval of observation are c%mpared in Table I, where ?, refers to true average invasion rates. Corresponding values obtained by the point-area method and segmentally continuous input functions are ri, pa and r,, BE, respectively. Enterohepatic Recirculation-A hypothetical drug is assumed to impose on the body the properties of the one-compartment open body model with a volume of distribution of 10 L, referenced to the drug in plasma. Plasma protein binding and red blood cell partitioning are assumed to be negligible, and the drug is not metabolized but excreted into bile and by glomerular filtration into urine. Renal clearance, referenced to plasma, and biliary clearance, referenced to blood (hematocrit: 0.5), are 120 and 260 mL/min, respectively, so that the hepatic extraction ratio is E = 0.2. The first-order rate constants are Iz, = 0.025, k,, = 0.012, and k 1 4 = 0.013 min-'. As discussed above, the drug is injected and bile is released 120 min after administration. During this interval, 49.4% of the initial dose accumulated in the bile, 80% of which was not subject to the first-pass effect. In the simulation, the rate constant of removal from compartment 0 was k,, = 0.05. The resulting plasma levels are shown in Fig. 2. They were analyzed by nonlinear regression assuming no drug input between open circles, which served to estimate the volume of distribution and the overall elimination rate constant. For the intervals between solid dots, segmentally continuous zero-order input was assumed and corresponding rates were determined. The relationship between average input rates within sampling intervals and the first-order input function used to generate the data is demonstrated in Fig. 3.
Discussion The model functions presented are a logical extension of classical linear pharmacokinetic models, which may, in fact, be treated as special cases. In the field of data analysis, they fill the gap between compartmental analysis using continuous input functions and numerical deconvolution methods. This does Journal of Pharmaceutical Sciences Vol. 74,No. 2,February 1985
1
139
h
$61
v
Ib
\
minutes Figure 2-Simulated time course of plasma concentrations with enterohepatic recirculation.
times related to gastric emptying and intestinal transport. If absorption phases are sufficiently long, the disappearance of drug from the absorption site may be accounted for by a sequence of first-order absorption processes, where the rate constants reflect, among other physiological variables, the permeability of the intestinal wall, which need not be constant over the full length of the gut. If data from an intravenous study are not available, experimental evidence may, in simple cases, be substituted to a certain extent by reasonable assumptions. Similar to the WagnerNelson method, an overall elimination rate constant can be estimated from data points in the terminal phase assuming that the one-compartment body model holds and that absorption becomes negligible early enough. Then, attributing an arbitrary value (e.g., 1)to the unknown volume of distribution, a staircase function of relative absorption rates may be computed for the early sampling intervals. If these rates are multiplied by a common factor so that the product with the total area under the staircase input is 100, the products of this factor with the fractional areas for the intervals correspond to the percentages of drug, referenced to the total amount absorbed, absorbed between pairs of observations. Besides the simulation and analysis of enterohepatic recirculation, the time course of plasma concentrations on repetitive or nonrepetitive multiple dosing at unequal intervals or with unequal single doses may be studied conveniently. Further applications may lie in the assessment of the reliability of parameter estimates, since the computation of sensitivity and dispersion matrices is simplified because partial derivatives are directly available, and in the dose optimization of sampling
schedule^.^^.^^
Appendix: Glossary minutes Figure 3-Simulated first-order input and reconstructed segmentally continuous zero-order input after enterohepatic recirculation. Key: (0) points used to estimate parameters of disposition model; (0)points used to estimate staircase input.
not imply that they can compete with the latter in terms of efficiency when used for curve fitting, since the time required to locate a minimum depends on the complexity of the problem, the goodness of the initial parameter estimates, the performance of the search algorithm, and the convergence criterion. The current unsatisfactory state may, however, be significantly improved by using more powerful iteration procedures, for which the availability of explicitly computed partial derivatives with respect to all model parameters is a promising feature. What makes the method attractive is its flexibility. The typical evaluation of absorption profiles proceeds in three steps: 1. First, the parameters of the disposition function are determined from plasma levels observed after intravenous administration of the drug. Invasion is treated as a zero-order process with known to, ri, and tl since the simplifying assumption of instantaneous input would require a more complicated treatment in this case. 2. Next, a staircase input function is obtained for data measured after extravascular administration, holding the parameters of the disposition function constant and equating sampling times with transition times. In this step, the point-area deconvolution method is a favorable alternative. 3. Statistically, staircase input functions are, however, not parsimonious with respect to the number of parameters estimated since one input rate is computed for each observation. Therefore, in a third step, the input function may be simplified systematically, replacing fixed transition times derived from the sampling schedule by physiologically meaningful transition 140
Journal of Pharmaceutical Sciences Vol. 74, No. 2, February 1985
C,
drug concentration in plasma dose D,.r apparent dose delivered to the intestinal absorption site after the lth emptying of the gallbladder sum of the first-order rate constants for removal of E, drug from the peripheral compartment j ; also, the apparent rate constant of return from the corresponding peripheral compartment of a mammillary model G ( t ) unit-impulse response function (synonyms: weighting function, system function) in,(t) input rate function for compartment 1 first-order rate constant of disapperance of the drug from the absorption site; apparent absorption rate constant apparent first-order absorption rate constant during the p t h absorption phase apparent first-order overall elimination rate constant referenced to the central compartment of a mammillary model first-order rate constant of drug transfer from compartment i to compartment j eigenvalues of a characteristic equation: observable hybrid rate constants rate of drug input, corresponding to constant infusion or other zero-order invasion process rate of drug input during the p t h invasion phase ith sampling time; in stepwise data evaluation: ith transition time time of the Ith emptying of the gallbladder after start of the experiment volume of distribution, one-compartment open body model volume of compartment 1,the central compartment of a mammillary model used for reference
D
xi
lp
amount of drug in compartment i amount of drug in compartment i after Ith emptying of gallbladder switch variable for truncation correction after the p t h absorption phase
References and Notes 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.
Teorell, T. Arch. Intern. Phurmucodyn. 1937,57, 205. Teorell, T. Arch. Intern. Pharmucodyn. 1 9 3 7 , 57, 226. Cutler, D. J. J. Phurmucokinet. Biophurm. 1978, 6, 227. Cutler, D. J. J. Pharmacokinet. Biopharm. 1978,6,243. Veng-Pedersen, P. J. Phurm. Sci. 1980,69, 298. Veng-Pedersen, P. J. Pharm. Sci. 1980,69,305. Veng-Pedersen, P. J. Phurm. Sci. 1980,69,312. Veng-Pedersen, P. J. Phurm. Sci. 1980,69,318. Wagner, J. G.; Nelson, E. J. Phurm. Sci. 1964,53,1392. Loo, J. C. K.; Riegelman, S. J. Phurm. Sci. 1 9 6 8 , 57,980. Vaughan, D. P.; Dennis, M. J. Phurmucokinet. Bwphurm. 1 9 8 0 ,
8. 83. 12. Vkghan, D. P.; Dennis, M. J. Pharm. Sci. 1978,67,663. 13. Chiou, W. L. J. Pharm. Sci. 1980,59,57. 14. Stepanek, R. “Praktische Analyse Linearer Systeme durch Fal-
tungsoperationen”; Geest u. Portig: Leipzig, 1976. 15. Benet, L. Z. J. Pharm. Sci. 1 9 7 2 , 6 1 , 532. 16. Vauahan, D. P.: Mallard, D. J. H.; Trainor, A.; Mitchard, M. Eur. J. C h i nPharmacol. . 1 9 7 5 , 8, 141. 17. Vaughan, D. P.; Trainor, A. J. Phurmacokinet. Bwpharm. 1975. 3 . 203. 18. Wolf, M.; Heinzel, G.; Koss, F. W.; Bozler, G. Arzneim.-Forsch. 1977,27(1), 900. 19. Heinzel, G.; Hammer, R.; Wolf, M.; Koss, F. W.; Bozler, G. Arzneim.-Forsch. 1977,27(1), 904. 20. Veng-Pedersen, P. J. Phurm. Sci. 1 9 7 8 , 6 7 , 187. 21. Sharney, L.; Wasserman, L. R.; Gevirtz, N. R. Am. J. Med. Electron. 1 9 6 4 , 4, 249. 22. Wagner, J. G. J.Phurmacokinet. Biophurm. 1976,4,395. 23. Vaughan, D. P.; Dennis, M. J. Pharmacokinet. Biophurm. 1 9 7 9 , 7. 511. 24. Suverkriip, R. J.Phurm. Sci. 1979,68,1395. 25. Nelder, J. A.; Mead, R. Comput. J. 1 9 6 5 , 7, 308. 26. Suverkriip, R. in “Pharmacokinetics During Drug Development: Data Analysis and Evaluation Techniques”; Bozler, G.; van Rossum, J. M., Eds.; G. Fischer: Stuttgart, 1982. 27. These programs are available from the author upon request. They are written in FORTRAN IV and are recorded on 1600 bpi tape. Prepayment of US $20.00 for tape, shipping, and documentation is required. - I
Journal of Pharmaceutical Sciences Vol. 74, No. 2, February 1985
141
partmental disposition functions can be used to analyze and simulate the time course of drug invasion into the central compartment of mammiliary models and certain physiologically relevant recirculating systems. The generalized model equations may be used to assess the reliability and intercorrelations of parameter estimates directly, since partial derivatives with respect to all model constants can be computed explicitly. In combination with curve-fitting algorithms, input functions identical to those of the point-area deconvolution can be obtained, although at the expense of significantly longer computing times. On the other hand, the range of potential applications goes far beyond the reconstruction of the time course of drug absorption. Smooth input functions are, in many cases, only a rough approximation of rather erratic absorption profiles; the time course of absorption is possibly subject to greater intra- and intersubject variability than most other processes of pharmacokinetic interest, but it is a prerequisite for the rationai design of dosage forms. The algorithms commonly employed to analyze drug invasion characteristics can be grouped into three classes: 1. Curve fitting of explicit convolutions of simple input functions with the disposition functions of compartmental models, which may be specified in the form of microscopic rate constants or as parameters of the general solution for plasma concentration profiles after impulse input, represent the first class. Typical examples are the functions given by Teorell,'~2 C ~ t l e r ,and ~ . ~Veng-Peder~en.~-' Generally suitable input functions should be flexible enough to describe the general characteristics of drug invasion in sufficient detail, yet so stiff that single measurements which do not fit into the overall pattern do not cause them to go astray. 2. Algorithms based on mass-balance considerations like the methods of Wagner and Nelson9 and Loo and Riegelman" are model independent with respect to input kinetics but rely on assumptions concerning the structure of the disposition model, since microscopic rate constants are used in the course of computations. Vaughan and Dennis" have demonstrated that the Loo-Riegelman method remains operative even if these assumptions do not strictly hold. 3. Numerical deconvolution procedures like the point-area method of Vaughan and Dennis," midpoint-impulse of Chiou,I3 and the trapezoidal approximation of Stepanek,I4 which is practicable even if a computer is not available, are derived from linear systems theory and do not imply any prerequisites with respect to the type of input kinetics nor specific properties of the model beyond its linearity and time independence. Experimental error is attributed directly to the input function and may give rise to temporarily negative apparent absorption rates in the vicinity of the peak, at the plateau, or on the downslope of the plasma level curve, but not while plasma levels rise markedly, unless error is large. 136
1 Journal of Pharmaceutical Sciences Vot. 74, No. 2, February 1985
The problem is aggravated if plasma levels after intravenous bolus injection are used directly as the weighting function. This is possible if equally spaced samples have been taken at the same times after intravenous administration and absorption. Under these conditions, most deconvolution algorithms assume a particularly simple and attractive form. The shortcut suggests itself in the case of the trapezoidal approximation, which is equivalent to the point-area method if fractional areas under the weighting function are computed by the trapezoidal rule. If this is done, the errors of two data sets are compounded and the outcome is uncertain. Sampling schedules employed in practice rarely coincide after intravenous administration and absorption, and samples are hardly ever taken at equal intervals. Therefore, in order to minimize the computing effort, data are frequently interpolated. At the same time, both the impulse response and the response function to be deconvoluted may be smoothed to make the input profile look more acceptable. This appears to be largely a cosmetic operation, which should be considered carefully, since systematic errors may be introduced even if the data are flawless.
Theoretical Section A general disposition function for linear mammillary models was given by Benet" and extended by Vaughan and cow o r k e r ~ ' ~to . ' ~include input into peripheral compartments. The matrix approach to Laplace transform solutions for linear multicompartment systems arbitrarily interconnected by firstorder transfer has been developed by Wolf, Heinzel, and coworker~'~.'~ and by Veng-Pedersen?' whose notation will be used below. Although, in the general case, input may occur at any time into any compartment with variable rates, only the particular situation where input is restricted to one compartment and enters the system as an impulse at time zero will be considered here, in order to define the weighting function of the system with respect to this compartment. The differential equations describing disposition in a dense n-compartment system with rate constants greater or equal to zero may be written:
+
n
c + ... +
+
j=l j#i
+ ...
+ kin
... -
n-1
c
j=1
0022-3549/85/0200-0136$0 1 .0010 0 1985, American PharmaceuticalAssociation
or briefly: x ' = (KT- Z)x
(2)
K Tis the n x n matrix with rate constants kGfor drug transfer from the ith to the j t h compartment, i = 1 , 2 , . . + ,n and j = 1, 2, . . . , n, where k, = 0 for i = j and Z is a diagonal matrix with elements Ej: n
E.= C k..
compartment model during and after infusion with rate ri, may be written:
In the case of first-order absorption, a lag time to may be taken into account, by which the continuous input response function is also shifted in time, and the first-step variable is defined as:
i=l i#j
With initial conditions xl.o = D and xi.0 = (i = 2, 3, . . .,n),the solution for x1 in the Laplace domain can be written:
D -fl(s)=
(s
+ Ej)
j-2 n
n (s +
In the definition of {I, a scale reduction due to the removal of drug from the absorption site in the time interval to+ has to be introduced, and we have:
(3) Ah)
k-1
where X k are the n eigenvalues of - Z. Dividing the inverse transform by the scaling factor V , , which relates amounts in this compartment to concentrations, the unit-impulse response function is:
Thus, the time course of drug concentrations in plasma for the one-compartment open body model with truncated first-order absorption is:
n
k=l
k#i
In pharmacokinetics, further restrictions with respect to the type of models are common: complex eigenvalues may arise if compartments are connected in loops, resulting in damped sinusoidal oscillations in open systems. Since this is hardly observed in reality, interest has been focused on mammillary models in the broad sense, where branched systems of peripheral compartments may be connected to the central compartment but lateral interconnections forming loops are absent.'l A solution of the form given above can always be derived starting from a model of given structure and known microscopic rate constants, but it is not possible to define both structure and values of rate constants unambiguously from the parameters E,,X i , and V , of the general solution. The problem of vanishing exponentials has been discussed by Wagner22and by Vaughan and who have shown that the number of models corresponding to any given system function is countably infinite if only plasma concentrations are measured. This family of models may, however, be mapped uniquely onto its simplest member, the n-compartment mammillary model with elimination only from the central compartment. Input Function-If a zero- or first-order drug input into a linear compartmental system exhibiting multiexponential input response is stopped abruptly, the time course of plasma concentrations can be represented as the sum of a continuous input response function and a truncation correction function of identical form.24In the case of zero-order infusion starting at time zero, the only differences are the sign of the corrective function and its displacement along the time axis by t l , the time when infusion stops. The notation becomes particularly Let 5; be simple if two-value step variables are intr~duced.'~ defined as: (5)
Then, the time course of plasma concentrations in the one-
If the invasion rate (rin)or the absorption rate constant (k.) do not drop to zero at time t l , but assume values of rin.2or k,.2,respectively, differing from the corresponding values of the first phase (which are now designated and k,.l), a second phase follows. The times at which rate constants are assumed to change their values may be called transition times. For zero-order invasion, a second step variable has to be defined, which changes its value at time t2:
Then, eq. 6 can be amended to read
Similarly, for first-order absorption:
For the case of m-phasic zero-order invasion, a general input function in the domain of the Laplace transform can be derived
The corresponding function for m-phasic first-order absorption is somewhat more complex, but still transparent:
x
[e-~ta-l
-
e-ka,(t.-t.-l)
e
-3t,
(14) l
The convolution product of these segmentally continuous Journal of Pharmaceutical Sciences Vol. 74. No. 2, February 1985
/ 137
input functions with m invasion phases and the general disposition function for n-compartment systems is readily farmed. After inverse transformation into the time domain, the plasma level function for multiphasic zero-order invasion reads:
i=l
- ti
with: (16)
Although the function looks unwieldy a t first glance, its structure is highly regular. The zero-order input response function consists of a time-independent steady-state term and n negative exponential terms. Both numerators and denominators of their coefficients are continued products formed after the same pattern. Note that the steady-state term may be converted to the same form by introducing a virtual constant Xo = 0. For each complete invasion phase, there is one continuous input and one corrective partial function. The corresponding plasma level function for segmentally continuous first-order input reads:
derivatives with respect to all model parameters can be computed explicitly. Enterohepatic Recirculation-The approach developed for mammillary models in the broad sense is also applicable to certain types of models with closed loops describing enterohepatic recirculation. Although the production of bile is a continuous process, bile flow in humans is not steady, as in some animal species. In humans, the emptying of the gallbladder is triggered by stimuli related to food intake. In the context of multicompartmental systems, the contents of the gallbladder may be considered external to the body, like the contents of the urinary bladder or of the gastrointestinal tract. In order to keep the following example simple, first-order transfers are assumed throughout (for the time being), with the only exception of bile release into the duodenum. Then, with one-compartment disposition kinetics and assuming the liver to be a compartment of vanishing volume, Fig. 1is a schematic representation of the model. As long as bile accumulates in the gallbladder, the switch between compartments 4 and 0 is open, and the following system of differential equations holds prior to the first emptying at time tb x; =
-(ko1
+ kz + ko4)Xo
x: =
kolX0
- (kl2 + k13 + k14)XI
x; =
kol xo
hzXi
xi =
ki3Xi
kQ4xo
x; =
(19)
If a dose D is introduced into compartment 1 as a bolus, we have initial conditions a t t = 0, x1 = D,xo = x2 = x3 = x4 o = 0. Immediately prior to the emptying at time tb the distribution of drug amounts in the system is:
-
n
x3
= ki3
x4 = k14
+E
kl4xl
1'"' 1"'
x1 dt
x1 dt
(20)
For the sake of simplicity, again it is assumed that the whole amount of bile accumulated up to t b . 1 is released in so short a time interval ( A t ) that the change of amounts of drugs in compartments 1,2,and 3 is negligible. Subsequently, the transport from compartment 4 to compartment 0 is inhibited and eq. 19 can be solved until the gallbladder is emptied the next e
(17)
1*1
with:
(18)
Because of their repetitive pattern of formation, these equations are readily implemented in the form of compact, yet structurally transparent, computer programs. Even their partial 138
Journal of Pharmaceutical Sciences Vol. 74, No. 2, February 1985
Figure 1-Compartmental model for enterohepatic recirculation. Key to compartments: (0) GI tract; (1) general circulation; (2) metaboiites; (3) drug excreted unchanged into urine; (4) gallbladder.
Table I-Comparison of Staircase Invasion Functions Obtained by Point-Area Deconvolution and Nonlinear Regression Using a Multiphasic Zero-Order Invasion Model with True Values
time at tb.2. Now, the initial conditions are:
0.0 1 .o
2.0 3.0 4.0 5.0 6.0 7.0 X4.tbl
(21)
=0
It is convenient to define the collective rate constants for drug elimination as k, = klz + kI3 + k14 and for its removal from the GI tract as k. = kl+ k2+ k4. The solution for x1 is obvious; superimposed upon the primary impulse response function: = D.e-%'
(22)
and with f = k,,,/k., we have for the amount of recirculated drug after the first release of bile:
0.0 4.25 5.31 5.06 4.37 3.60 2.89 2.30
0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0078
-
0.4912 0.2463 0.1230 0.0622 0.0311 0.0147 0.0077
-
0.4912 0.2463 0.1230 0.0622 0.0306 0.0155 0.0075
-
computer-generated data with and without random error added and from experimental data, taking the times of observations as transition times. The staircase input functions arrived at were virtually identical with those obtained by the point-area deconvolution method of Vaughan and Dennis. The unit-impulse response function given in their paper:"
G ( t ) = 5.Oe-O.*' + 5.0e-'.'
(24)
can be formulated according to eq. 4 as:
Its convolution with the first-order input function: inl(t) = 0.693e-0693' The total amount of drug in compartment 1 is XI = x1 , + x1 r ~ . From time t b on, the intermediate reservoir in compartment 4 is fed from two sources: besides the biliary accumulation of drug cleared from the general circulation, a first-pass process may contribute significantly. Looking only at compartment 1 (drug in plasma) and possibly at compartment 3 (unchanged drug excreted into urine), the system behaves as if a second dose D, has been introduced into compartment 0 at time tb 1. The extension to more than one instance of bile release and the multicompartmetal models of drug disposition is conceptually simple. Several simplifications have been introduced into the example for the sake of clarity, which may be dropped in the general case. Since compartment 0 is considered to be external and a predecessor of compartment 1, the linearity of the model is maintained, as far as distribution and elimination are concerned, if first-order processes of drug absorption, presystemic metabolism, and presystemic biliary secretion are replaced by arbitrary transfer functions, all greater than or equal to zero at any time. Moreover, instantaneous and complete emptying of the gallbladder is not a strict requirement if interest is focused on the time course of absorption of recirculated drug into compartment 1. In order to define the disposition characteristics of a multicompartment system adequately after intravenous administration, i.e., to estimate the parameters E, and h, of its general solution for its circulation compartment with reasonable precision, there must only be sufficiently long periods in the initial and terminal phases during which release of bile and reabsorption are negligible.
Experimental Section For the purpose of data analysis, a subroutine incorporating eqs. 15-18 has been combined with the Nelder-Mead simplex algorithm for minimization of the sum of (weighted) residual squares.26This method, although less efficient than other optimization procedures, has been chosen because it allows the user to fix the values of known model constants conveniently by specifying zero initial variation steps for these parameters. Comparison wit h Point-Area Deconvolution-Segmentally continuous zero-order input rates were estimated from
(26)
was analyzed both by deconvolution and curve fitting, holding the values of V ] , EP,XI, and hz constant. The zero-order invasion rates obtained. for each interval of observation are c%mpared in Table I, where ?, refers to true average invasion rates. Corresponding values obtained by the point-area method and segmentally continuous input functions are ri, pa and r,, BE, respectively. Enterohepatic Recirculation-A hypothetical drug is assumed to impose on the body the properties of the one-compartment open body model with a volume of distribution of 10 L, referenced to the drug in plasma. Plasma protein binding and red blood cell partitioning are assumed to be negligible, and the drug is not metabolized but excreted into bile and by glomerular filtration into urine. Renal clearance, referenced to plasma, and biliary clearance, referenced to blood (hematocrit: 0.5), are 120 and 260 mL/min, respectively, so that the hepatic extraction ratio is E = 0.2. The first-order rate constants are Iz, = 0.025, k,, = 0.012, and k 1 4 = 0.013 min-'. As discussed above, the drug is injected and bile is released 120 min after administration. During this interval, 49.4% of the initial dose accumulated in the bile, 80% of which was not subject to the first-pass effect. In the simulation, the rate constant of removal from compartment 0 was k,, = 0.05. The resulting plasma levels are shown in Fig. 2. They were analyzed by nonlinear regression assuming no drug input between open circles, which served to estimate the volume of distribution and the overall elimination rate constant. For the intervals between solid dots, segmentally continuous zero-order input was assumed and corresponding rates were determined. The relationship between average input rates within sampling intervals and the first-order input function used to generate the data is demonstrated in Fig. 3.
Discussion The model functions presented are a logical extension of classical linear pharmacokinetic models, which may, in fact, be treated as special cases. In the field of data analysis, they fill the gap between compartmental analysis using continuous input functions and numerical deconvolution methods. This does Journal of Pharmaceutical Sciences Vol. 74,No. 2,February 1985
1
139
h
$61
v
Ib
\
minutes Figure 2-Simulated time course of plasma concentrations with enterohepatic recirculation.
times related to gastric emptying and intestinal transport. If absorption phases are sufficiently long, the disappearance of drug from the absorption site may be accounted for by a sequence of first-order absorption processes, where the rate constants reflect, among other physiological variables, the permeability of the intestinal wall, which need not be constant over the full length of the gut. If data from an intravenous study are not available, experimental evidence may, in simple cases, be substituted to a certain extent by reasonable assumptions. Similar to the WagnerNelson method, an overall elimination rate constant can be estimated from data points in the terminal phase assuming that the one-compartment body model holds and that absorption becomes negligible early enough. Then, attributing an arbitrary value (e.g., 1)to the unknown volume of distribution, a staircase function of relative absorption rates may be computed for the early sampling intervals. If these rates are multiplied by a common factor so that the product with the total area under the staircase input is 100, the products of this factor with the fractional areas for the intervals correspond to the percentages of drug, referenced to the total amount absorbed, absorbed between pairs of observations. Besides the simulation and analysis of enterohepatic recirculation, the time course of plasma concentrations on repetitive or nonrepetitive multiple dosing at unequal intervals or with unequal single doses may be studied conveniently. Further applications may lie in the assessment of the reliability of parameter estimates, since the computation of sensitivity and dispersion matrices is simplified because partial derivatives are directly available, and in the dose optimization of sampling
schedule^.^^.^^
Appendix: Glossary minutes Figure 3-Simulated first-order input and reconstructed segmentally continuous zero-order input after enterohepatic recirculation. Key: (0) points used to estimate parameters of disposition model; (0)points used to estimate staircase input.
not imply that they can compete with the latter in terms of efficiency when used for curve fitting, since the time required to locate a minimum depends on the complexity of the problem, the goodness of the initial parameter estimates, the performance of the search algorithm, and the convergence criterion. The current unsatisfactory state may, however, be significantly improved by using more powerful iteration procedures, for which the availability of explicitly computed partial derivatives with respect to all model parameters is a promising feature. What makes the method attractive is its flexibility. The typical evaluation of absorption profiles proceeds in three steps: 1. First, the parameters of the disposition function are determined from plasma levels observed after intravenous administration of the drug. Invasion is treated as a zero-order process with known to, ri, and tl since the simplifying assumption of instantaneous input would require a more complicated treatment in this case. 2. Next, a staircase input function is obtained for data measured after extravascular administration, holding the parameters of the disposition function constant and equating sampling times with transition times. In this step, the point-area deconvolution method is a favorable alternative. 3. Statistically, staircase input functions are, however, not parsimonious with respect to the number of parameters estimated since one input rate is computed for each observation. Therefore, in a third step, the input function may be simplified systematically, replacing fixed transition times derived from the sampling schedule by physiologically meaningful transition 140
Journal of Pharmaceutical Sciences Vol. 74, No. 2, February 1985
C,
drug concentration in plasma dose D,.r apparent dose delivered to the intestinal absorption site after the lth emptying of the gallbladder sum of the first-order rate constants for removal of E, drug from the peripheral compartment j ; also, the apparent rate constant of return from the corresponding peripheral compartment of a mammillary model G ( t ) unit-impulse response function (synonyms: weighting function, system function) in,(t) input rate function for compartment 1 first-order rate constant of disapperance of the drug from the absorption site; apparent absorption rate constant apparent first-order absorption rate constant during the p t h absorption phase apparent first-order overall elimination rate constant referenced to the central compartment of a mammillary model first-order rate constant of drug transfer from compartment i to compartment j eigenvalues of a characteristic equation: observable hybrid rate constants rate of drug input, corresponding to constant infusion or other zero-order invasion process rate of drug input during the p t h invasion phase ith sampling time; in stepwise data evaluation: ith transition time time of the Ith emptying of the gallbladder after start of the experiment volume of distribution, one-compartment open body model volume of compartment 1,the central compartment of a mammillary model used for reference
D
xi
lp
amount of drug in compartment i amount of drug in compartment i after Ith emptying of gallbladder switch variable for truncation correction after the p t h absorption phase
References and Notes 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.
Teorell, T. Arch. Intern. Phurmucodyn. 1937,57, 205. Teorell, T. Arch. Intern. Pharmucodyn. 1 9 3 7 , 57, 226. Cutler, D. J. J. Phurmucokinet. Biophurm. 1978, 6, 227. Cutler, D. J. J. Pharmacokinet. Biopharm. 1978,6,243. Veng-Pedersen, P. J. Phurm. Sci. 1980,69, 298. Veng-Pedersen, P. J. Pharm. Sci. 1980,69,305. Veng-Pedersen, P. J. Phurm. Sci. 1980,69,312. Veng-Pedersen, P. J. Phurm. Sci. 1980,69,318. Wagner, J. G.; Nelson, E. J. Phurm. Sci. 1964,53,1392. Loo, J. C. K.; Riegelman, S. J. Phurm. Sci. 1 9 6 8 , 57,980. Vaughan, D. P.; Dennis, M. J. Phurmucokinet. Bwphurm. 1 9 8 0 ,
8. 83. 12. Vkghan, D. P.; Dennis, M. J. Pharm. Sci. 1978,67,663. 13. Chiou, W. L. J. Pharm. Sci. 1980,59,57. 14. Stepanek, R. “Praktische Analyse Linearer Systeme durch Fal-
tungsoperationen”; Geest u. Portig: Leipzig, 1976. 15. Benet, L. Z. J. Pharm. Sci. 1 9 7 2 , 6 1 , 532. 16. Vauahan, D. P.: Mallard, D. J. H.; Trainor, A.; Mitchard, M. Eur. J. C h i nPharmacol. . 1 9 7 5 , 8, 141. 17. Vaughan, D. P.; Trainor, A. J. Phurmacokinet. Bwpharm. 1975. 3 . 203. 18. Wolf, M.; Heinzel, G.; Koss, F. W.; Bozler, G. Arzneim.-Forsch. 1977,27(1), 900. 19. Heinzel, G.; Hammer, R.; Wolf, M.; Koss, F. W.; Bozler, G. Arzneim.-Forsch. 1977,27(1), 904. 20. Veng-Pedersen, P. J. Phurm. Sci. 1 9 7 8 , 6 7 , 187. 21. Sharney, L.; Wasserman, L. R.; Gevirtz, N. R. Am. J. Med. Electron. 1 9 6 4 , 4, 249. 22. Wagner, J. G. J.Phurmacokinet. Biophurm. 1976,4,395. 23. Vaughan, D. P.; Dennis, M. J. Pharmacokinet. Biophurm. 1 9 7 9 , 7. 511. 24. Suverkriip, R. J.Phurm. Sci. 1979,68,1395. 25. Nelder, J. A.; Mead, R. Comput. J. 1 9 6 5 , 7, 308. 26. Suverkriip, R. in “Pharmacokinetics During Drug Development: Data Analysis and Evaluation Techniques”; Bozler, G.; van Rossum, J. M., Eds.; G. Fischer: Stuttgart, 1982. 27. These programs are available from the author upon request. They are written in FORTRAN IV and are recorded on 1600 bpi tape. Prepayment of US $20.00 for tape, shipping, and documentation is required. - I
Journal of Pharmaceutical Sciences Vol. 74, No. 2, February 1985
141