Time Varying Apparent Volume of Distribution and Drug Half-Lives Following Intravenous Bolus Injections

We present a model that generalizes the apparent volume of distribution and half-life as functions of time following intravenous bolus injection. This generalized model defines a time varying apparent volume of drug distribution. The half-lives of drug remaining in the body vary in time and become longer as time elapses, eventually converging to the terminal half-life. Two example fit models were substituted into the general model: biexponential models from the least relative concentration error, and gamma variate models using adaptive regularization for least relative error of clearance. Using adult population parameters from 41 studies of the renal glomerular filtration marker 169Yb-DTPA, simulations of extracellular fluid volumes of 5, 10, 15 and 20 litres and plasma clearances of 40 and 100 ml/min were obtained. Of these models, the adaptively obtained gamma variate models had longer times to 95% of terminal volume and longer half-lives.


Introduction
The apparent volume of distribution (V d ) is an important pharmacokinetic parameter that relates drug plasma concentrations to the amount of drug in the body and is important for drug loading dose and maintenance dose calculations [1,2,3]. Following an intravenous bolus of drug, the volume of distribution of a drug varies with time. This redistribution occurs in two distinct phases with very different time scales, a vascular phase and a washout or dilution phase. During the vascular phase, a drug mixes throughout the blood's plasma volume with a time scale of seconds or minutes. During the dilution phase, hydrophilic drugs distribute into the body's interstitial fluids with a time scale of hours or days. While most drugs are actively redistributing from plasma into the body tissues, the decrease of plasma concentrations is largely due to that redistribution as opposed to actual drug elimination [2]. In the following, we assume first order kinetics, i.e., that drug elimination is proportional to its concentration, which is the most common drug kinetic. The most commonly calculated volumes of distributions are the apparent volume of drug distribution immediately after bolus intravenous injection, i.e., at time zero, (V 0 ), the terminal apparent volume of distribution following bolus intravenous administration (V area ) and the expected volume of distribution, called V E here, and often called Vss in the literature. V SS for a constant infusion experiment is the terminal apparent volume of distribution, analogous to V area in a bolus model [2,3]. For bolus experiments, in place of V SS we use the term V E for the expectation of a physical volume of distribution of drug, and unlike V SS , V E is invariant between constant infusion and bolus experiments.
Time varying apparent volume of distribution models based on exponential washout models for plasma concentration, have been developed previously [2,4,5]. Niazi [4] developed a temporally variable, apparent volume of distribution model for sum of exponential (SET) functions based on conservation of mass. This variable volume model implies an explicit relationship between redistribution and the rate of volume of drug distribution expansion in time. Such a model specifies when the drug volume reaches a particular size relative to its terminal apparent volume of distribution. In so doing, such a variable volume model could potentially provide unanticipated new information about time based tissue drug effects for tissue metabolism and/or eliminations. This type of model also has implications for effects of drugs on body tissues, i.e., therapeutic, toxic, or radiation exposure effects [3].
A variable volume model could be used to calculate the optimal instantaneous dose that will produce a desired concentration at the effector site without an overload of the drug. Here we present a general time-dependent apparent volume of distribution model, based on mass conservation, which could be used with almost any concentration washout fit function. To investigate how the implications of the general model relate to specific fit function models, we substitute biexponentials (E2) and gamma variates (GV) [6,7,8] into the general washout model.

Theory
The variable volume of distribution Plasma clearance (CL) is a drug's elimination rate, M 0 (t), divided by its corresponding plasma concentration, C(t), at any time t [9], where M(t) is total drug mass in the body at time t and the prime indicates differentiation. Rearrangement of Eq (1) gives the drug's elimination rate, Integrating both side of Eq (2) from time 0 to t gives mass remaining in the body at time t to be where M(t!1) = 0.When t = 0 and M(t = 0) = D, the dose, from equality of the left and right hand functions of Eq (3), one solves for CL to obtain the well-known area under the curve, AUC, definition of CL, Conservation of mass for a concentration in an apparent (homogenous plasma concentration) volume of distribution, V d (t), at time t, is given by The typical pharmacokinetic model specifies the washout of an impulse-response. The impulse is the entire initial dose M(0) = D distributed in some initial volume no matter how small. M(t) is monotonic non-increasing (mass conservation) and V d (t) monotonic increasing. Thus, C(t) is monotonic decreasing, and has a maximum value, no matter how large, at time is zero. Substitution of Eq (5) into Eq (3) yields the volume of distribution at any time t as Expected volume of distribution Expected volume of distribution (V E ) is the expected value of physiological volume of distribution of drug in the plasma and tissues of a concentration model. V E inherits the statistical properties for expectation from its distribution density function model, p(t), For example, if the distribution density function, p(t), has a mean value expectation, E½T ¼ t, we can then calculate the expected volume of drug distribution from that mean, which for a time series is called a mean residence time (MRT ¼ t) [10,11]. MRT is the average value or first moment of a time-series density function.
where P is the cumulative density function of p. To be useful, Eq (8) assumes that drug elimination rate is first order, i.e., that elimination is not dose dependent [12]. If a density function has no moments, MRT is undefined but there may still be a location that characterizes the data, for example, a Cauchy distribution has a stable median [T]. Indeed, for a distribution the median can be a better measure of location than the expected value, e.g., for some values of the beta distribution. For some density functions, E[T], is indeterminate for heavier tails, e.g., for the Pareto distribution.
We suggest that E[T] for a constant infusion experiment is also from p(t), and V E is from Eq (7). However, p(t) is the impulse-response to an instantaneous bolus, i.e., to a Dirac delta, and to obtain a p(t), we account for the integration of constant infusion as PðtÞ ¼ ð t 0 pðtÞdt, the cumulative density of p(t). We then suggest fitting D R CL PðtÞ ¼ C SS PðtÞ to constant infusion data, where D R is the dose rate of constant infusion, D R /CL is the terminal concentration of the infusion experiment usually called C SS . Once we have a p(t), its derivative is p(t).
Half-life of the drug, t 1/2 The half-life of any drug's mass in the body at time t can be written as where for a marker cleared only by the kidney this half-life corresponds to the urinary clearance half-life. Substituting Eqs (2 & 5) into Eq (9) yields an expression for half-life of a variable apparent volume However, the urinary clearance of kidney marker half-life is not, as a function of time, the same as the half-life of plasma concentration. To find the latter, we first substitute Eq (1) into the derivative of Eq (5) and rearrange terms to obtain which allow the half-life of plasma concentration as a function of time to be written as Given that V d (t!1) = V area (apparent volume at terminal phase), and V d 0 (t!1) = 0, a unique terminal half-life, t 1/2 (1), is the limiting value for both Eqs (10) and (12), i.e., for mass (e.g., renal) clearance and plasma concentration half-lives where at all other times, comparing Eqs (10 & 12), t 1/2 −M(t) > t 1/2 −C(t), such that before the terminal half-life becomes established, plasma concentration dilutes faster than mass clearance. To show this, Eqs (2, 5 & 11) are combined and terms rearranged to yield an explicit relationship between plasma dilution rate and mass clearance rate, such that the relative rate of plasma dilution is the relative rate of mass clearance plus redistribution, where redistribution is the relative rate of volume dilution. To see the effect sizes and their timing, we next substitute specific fit functions into the general model equations above, and later apply them to subject data.

Exponential solutions to the variable volume and half-life equations
In the exponential model, plasma concentration is given by sum of exponentials, where c i and λ i are the coefficients of i th exponential term. Substituting the sum definition of C exp (t) from Eq (15) into Eq (3) yields a sum of exponential terms solution to drug mass remaining in the body at time t, Similarly, substitution of the C exp (t) sum, Eq (15) into the general form for V d (t),Eq (6), yields, Substitution of t = 0,1 into this equation allows us to specify the initial (V 0 ) and final apparent volumes of distribution V area , In a compartmental model from sums of exponential terms, V 0 is identical to the volume of the central compartment V c .
Substitution of Eq (15) into Eq (4) yields the CL for SET model, Substitution of Eq (15) into Eqs (7 & 8) and performing the integration yields a steady-state volume of distribution (V SS ) for SET models, which in our more general context, we call V E , Substituting of the SET V d (t), Eq (17), into general half-life Eq (10) yields the mass clearance half-life for SET functions at time t, The mass clearance half-life would equal the renal clearance half-life for a GFR marker.

A gamma variate solution to the variable volume and half-life equations
Regularized GV functions are of interest because they have been previously shown to require one-half the sampling time (4 h) needed for numerical integration (8 h) to obtain precise and accurate CL-values in a large retrospective series [6]. The plasma concentration as a function of time can be modelled by gamma variate (GV) function, where α, β and K are the three parameters of a GV function. Note that α < 1 is not a constraint, and there is so far only one published method of consistently obtaining α < 1 values without using constraints [6,7,8,13]. Substitution of C GV (t) Eq (23) into general Eq (3) yields the gamma variate solution to drug mass remaining in the body at time t, where Γ(α) and Γ(α, βt) are respectively, the gamma function and the upper incomplete gamma function. Substitution of C GV (t) Eq (23) into Eq (6) and performing the indicated integration yields This equation is monotonically increasing only when α < 1, which is not a constraint for obtaining α-values, rather the meaning is that the upper limit for physicality of α-values is one. In that case, the initial (V 0 ) and final (V area ) apparent volume of distribution in the limit as time goes to 0 and to infinity are respectively given by Substitution of Eq (23), the C GV (t), into general Eq (4) and performing the indicated integration yields an expression for CL of a GV model Substitution of Eq (23) into Eqs (8 & 7) yields the physical volume of distribution for GV model, Note that from this and Eq (27) Gða; btÞ 1 À Gða;btÞ This limit exists except for α = 1, while β > 0. This latter does not occur for the GV solutions used here, which yield β!0+ when α!1−, i.e., concentration flatlines rather than become exponential [6]. When 0 < α < 1, V d is monotonic increasing and denominators of Eq (30) are positive valued. insures that M(t) < ε for V d (t ε < t < T), where T is a sufficiently large but finite time. Thus, V d (t) becomes as mass depleted as desired for t sufficiently large but less than some T-large, i.e., before V d converges to V area at infinite time. Consequently, the single bolus experiment GV model's drug depleted physical volume is written V E . That is, the GV washout model steady state is trivial and empty. The restriction on Eq (30) is that the GV not be an exponential. Readers who use hazard rates to classify tail heaviness of distributions may find this confusing. Hazard rate classification of tail heaviness is inexact and actual terminal tail areas compare as survival functions. From survival function ratios, gamma distributions with α > 1 have lighter than exponential tails, and for α < 1, i.e., the general case here, gamma distributions tails are heavier than exponential.
Substituting V d (t) from Eq (23) into Eq (10) yields the half-life of the drug mass remaining in the body Methods To test the time varying volume of distribution model, concentrations verses time curves for four V E values (10, 15, 20 and 25 L) at CL of 40 and 100 ml/min were simulated for biexponential (E2) and the gamma variate (GV) models. E2 and GV parameters were computed by using prior published data as follows. Data were used here from a prior study of 41 plasma concentration time samplings following intravenous bolus 169 Yb-DTPA (ytterbium diethylenetriaminepentaacetate) [14]. In this population, patients were given an antecubital IV bolus injection of 1.85 MBq of 169 Yb-DTPA. Eight blood samples were taken at 10, 20, 30, 45, 60, 120, 180 and 240 min after injection. Plasma clearance CL and V E for E2 and GV functions were calculated for all 41 patients. E2 functions were random search fit using relative-concentration weighted regression that preferentially fits longer times giving better results than from the use of ordinary least squares [15,16]. The Tk-GV algorithm uses Tikhonov regularization to increase covariance between the GV function and the time-samples by minimizing the AUC error over the entire interval from t = 0 to 1. This minimizes the relative error of plasma clearance [6]. The GV functions' three parameters were obtained from the Tk-GV method. An important point is that the Tk-GV method uses adaptive smoothing and without this feature the resulting PK parameter GV model results will be erratic [6,8]. A Windows compatible Tk-GV software application is available to individual researchers (i.e., not institutions) free of charge from the corresponding author.
Parameters for desired V E and CL values were obtained by interpolating parameters obtained from above fitted curves. Using computed parameters, volume of distribution, drug mass remaining in the body and drug half-lives as a function of time were plotted for four V Evalues at CL of 100 and 40 ml/min. We have ignored some of the different morphological implications of the different modelling in order to generate comparison curves. That is, CL and V E were assigned the same specific values, e.g., 100 ml/min and 20 litres, respectively, for both the E2 and Tk-GV models. For actual cases, that will not occur within the same hypothetical subject, and the same case between model CL-and V E -values differ slightly. A statistical analysis of the parameters for the reference population used appears in the Result Section below.  (Fig 2a and 2d) is increased with the increase of volume of distribution. On the other hand, Fig 2b and 2e show that dose in the body reduces at a faster rate for higher CL. In both Fig 1c and 1f, half-life is increased with the increase of the fluid volume. Furthermore, comparing same two figures (Fig 1c and 1f) at each V E shows that decrease of clearance tends to increase of the half-life.  When we consider the time to achieve 95% of V area , increase of V E increases the time to achieve 95% of V area . Fig 3b and 3e show that decreased dose remaining in the body occurred sooner for the larger CL-values, while smaller CL-value patients took more time to clear drugs from the body with an additional delay to achieve a given mass excretion seen with the increase of V E . Fig 3c and 3f show that half-lives were longest for the smaller CL-values. Furthermore, half-life increased for increasing V E . For comparison of the E2 and GV methods, Table 1 lists the time for V d to reach 95% of V area and drugs terminal half-life following intravenous boluses. Both time values shown in the Table 1 increase with the increase of the V E as well as with the decrease of the CL. Even though we selected equal CL and V E parameter values for our demonstration E2 and GV models, the GV models still predicted longer times for V d to reach 95% of V area as well as longer terminal half-lives.

Biexponential (E2) model
For same-case concentration data, it has been shown that the E2 and GV (Tk-GV) models will predict significantly, if slightly different CL and V E values [6,7,8]. To see the errors for these differences for the subject cases herein, bootstrap was used to avoid exaggerating the E2 model 95% confidence interval (CI) and standard deviation (SD) results that would occur if one spuriously assumed normal distribution conditions especially for the early exponential function's population parameters, which latter were found to be not normally distributed. For the Tk-GV method, the resulting gamma variate parameters were quasi-normally distributed, which allowed for the routine Tk-GV program SD-values to be used for calculating the CIs from their more normal distribution properties. Table 2 shows that the CIs are larger for E2 than for Tk-GV.
Some of the (many) parameters that were significantly different between models are summarized in Table 3. Note that V E , V area , t 1/2 and MRT are greater for Tk-GV than for E2, and that CL and CV(V E ) are lesser and as well β < λ 2 . Thus, although for easy comparison of families of curves in the figures, CL and V E , were assigned the same specific values for both the E2 and Tk-GV models, we were, in effect, contrasting two hypothetical subjects with slightly different body morphologies.
Fit quality for the Tk-GV and E2 models has been presented elsewhere for this adult data setsee Fig 1 and surrounding text of [6]. Herein, normalized covariances, i.e., correlations, between the logarithms of the time-samples and of each model evaluated at those times were considered indexed to the fidelity of trending between the data and the two models. These correlations had a median r-value for Tk-GV of 0.99840 and of 0.99832 for E2, with a two-tailed Wilcoxon p = 0.27 suggesting no significant difference in the overall quality of how each model trended with the data. However, even comparing r-values is somewhat misleading. The Tk-GV algorithm smooths its model curve to stabilize AUC on the 0 to infinite time interval, which time interval is unrelated to calculation of its r-value with the data, whereas the weighted-E2 models were fitted to concentration in a fashion numerically identical to the time-sample treatment used to compute r-values, such that the lack of significant improvement of r-values for weighted-E2 versus Tk-GV fitting does not increase our confidence in E2 model behaviour; weighted-E2 modelling was given an unfair advantage to outperform Tk-GV modelling, and did not do so.

Discussion
In both the biexponential (E2) and gamma variate (GV) variable volume models, smaller values of V E approached 95% of their terminal V d -values sooner than those with relatively greater  V E values, see Table 1. Figs 2b and 3b show that an increased fluid volume was associated with prolonged drug-levels at late time. These results are generally consistent with the known delays in redistribution from excess body fluid and in advanced renal insufficiency [7,17,18]. Furthermore, both the E2 and GV models demonstrated that drug body burden decreased sooner for show that for both model types, longer half-lives occurred for larger V E value patients in both CL of 100 and 40 ml/min. Drug elimination depends on the amount of drug delivered to the organs where excretion occurs. Therefore, the fraction of the drug residing in the plasma is a determinant of drug elimination. When V E is large, much of the drug is in extraplasmic space and is largely unavailable to feed any excretory organs. This is consistent with the herein observed increased V E values associated with increased drug half-lives such that drug remained in the body for longer times. For all four V E values at CL of 100 and 40 ml/min of the families of curves shown, the time to achieve a high percentage of the terminal apparent volume of distribution and terminal halflife were longer for the GV models than for the E2 models (Table 1). In this study, the E2 models achieved their terminal apparent volumes of distribution as well as terminal half-lives of drug mass within 2 hours for CL of 100 ml/min and within 16 hours for CL of 40 ml/min. On the other hand, GV model had not been reached the terminal half-life of mass even within 24 hours following the intravenous bolus of drug for all CL and V E values. The results in the Table 1 also show that GV model predicted longer terminal half-lives than the terminal halflives from the E2 model. A sum of exponential terms model reduces to a monoexponential for late time, and monoexponentials inherit their rate constants from the local data region fit while GV model extrapolate a terminal function by inheriting only a shape from a local data region.
Time samples of unmixed blood before arrival of a bolus in the sampling site are not used for washout modelling. Rather, the concentration fit function's area under the curve of a washout model is obtained by back-extrapolation from peripheral vein samples obtained when vascular mixing is already advanced. The variable volume approach used here allows any number of summed exponential terms to be directly substituted into the single variable volume model without invoking compartments, and, the compartmental approach is merely another explanation for the same SET fit equations. The physical ambiguity for SET models is nothing new, for example, Niazi admixed mammalian compartmental and variable volume models without distinguishing between them and Zierler identified 10 different physical configurations for any E2 model [4,19]. Using a variable volume model, an apparent volume of uniform drug concentration that increases in time becomes possible. That is, the variable drug volume represents the apparent volume, containing drug at any particular time, which contrasts starkly with the corresponding compartmental concept of a constant V area , that only corresponds to the apparent volume of distribution at infinite time in a variable apparent volume of distribution model. In washout modelling, E2 model usage implies instant mixing of drug in the central compartment, generally composed of plasma and additional spaces into which the drug distributes extremely rapidly and therefore, E2 models predict a finite apparent volume of distribution at time zero. Subsequently, the drug distributes from the central compartment to the other body spaces, e.g., into a peripheral compartment, and drug elimination commences simultaneously with the drug administration. In E2 models, the distribution process completes after a few hours [20]. For a two-compartment model, there is always a time at which V d transiently coincides with the value Vss and the further decrease of the plasma concentration of the drug is usually assumed to be only from drug elimination. The top of Fig 4 shows a schematic diagram of E2 model of drug distribution following bolus injection of a drug.
For GV models, the initial volumes of distribution are zero from Eq (26). The assumption of vanishingly small actual volume for the drug at the time of drug administration is realistic, V d is approximately zero at time zero, after which the drug distributes throughout the body. Other evidence to this effect is that back-extrapolated early concentrations of 51 Cr-EDTA, a GFR marker, have been observed to be more frequently of a linear logarithm of time shape, which functions set C(0) = 1, than of monoexponential shape, which do not [21]. Bottom Fig 4 shows a schematic diagram of GV model of drug distribution after bolus injection of a drug.
In prior work, in a 412 case series with both children and adults, some of us showed more precise and accurate CL-values for Tk-GV versus the less accurate E2 model compared to 8 h numerical integration controls for the slightly different GFR marker 99m Tc-DTPA [8]. In particular, the 4 h four sample E2 CL-values were significantly inflated by 4.9% compared to controls, and 4 h four sample Tk-GV CL-values were 0.5% (insignificantly) greater than eight sample, 8 h controls, for a net 4.2% CL difference between the 4 h E2 and Tk-GV experiments. In the current series, the 4 h, eight sample E2 CL-values were (significantly) 3.8% greater than the 4-h eight sample Tk-GV CL-values. The E2 model could also be drawn as a variable volume model in which case a scale factor α exp = V E /V d (1) < 1 would define the physical volume at time t to be α exp V d (t). Similarly, for the variable volume adaptively obtained GV model, one can define α = V E /V d (1) < 1, and an expanding physical volume αV d (t). Note, both α exp and α are constants at all times for their respective models. The term V SS can be confusing because 1) V SS implies that V E is always a steady state volume, which is not the case as the GV model αV d (t) < V E is concentration depleted at late time, see Eq (30). 2) V SS implies that V E only exists at t = 1, whereas V E is defined all of the time, i.e., on t = [0,1) by Eqs (8 & 7). Finally, 3) V SS implies an expected physical volume of distribution for sums of exponential term bolus models, and the apparent volume of distribution for a constant infusion experiment, whereas V E applies to more models as the expected volume of physical distribution of a drug for both the bolus and constant infusion experiments. In this series, the expected volumes of distribution results were significantly less variablesee Tables 2 and 3-for Tk-GV than for E2. Indeed, for all the tested quantities, there were no result groupings for the Tk-GV parameters that were more variable than those for the E2 models, although there are many possible comparisons and admittedly not all possible parameter comparisons were calculated, and some tested differences were not highly significantly different. Thus, it would appear that in so far as we were able to test that assertion, the Tk-GV model PK parameter results were more precise than the E2 model. Given this finding and prior better accuracy results [8], the 5 to 9 times fold longer 95% of V area times seen for the Tk-GV models are plausible. This may have important implications for future pharmacokinetic work.

Conclusion
We have presented initial examples of the power of the general approach to variable volume modelling using a general time dependent volume of distribution model based on mass conservation that sets time varying drug body mass as proportional to volume of distribution. Terminal half-life depends on both expected drug distribution volume and the clearance. Furthermore, two commonly used model functions were used to illustrate behaviour of the general model. Exploration of other yet-undescribed fit function models and the modelling of drugs having more complicated pharmacokinetics is left for future work.