^{1–6}Although a physiologically based pharmacokinetic model may be considered an accurate description of drug disposition from the start of drug administration, difficulty in collecting the extensive data for the model development and the fact that such models are not readily fitted to clinical pharmacokinetic-driven drug infusion caused by their complexity, have limited their applicability into pharmacokinetic studies.

^{5,7}More experimental approaches have been described to reveal the delay between drug administration and appearance of the drug at the site of arterial blood sampling. Standard lag time models (LAG model), which represent the time shift of dosing as if the drug was in fact administered at a later time, and analytical solutions of the TRANSIT compartment model (TRANSIT model), which depicts a multiple-step process represented by a chain of presystemic compartments and differs from the LAG model by the presence of one additional parameter, a transit rate constant, have been used by various authors.

^{8}Others have used more sophisticated recirculatory models to analyze arterial drug concentrations

*versus*time using estimated cardiac output from the first pass area under the curve in animal research.

^{7}

_{e0}, characterizes the drug transfer between the plasma and effect site.

^{9}When using a combined pharmacokinetic-dynamic model to describe the complete dose-response relationship, misspecification in the pharmacokinetic model over the first few minutes will result in a less accurate interpretation of the plasma effect-site equilibration.

^{6}

^{6,10}

#### Materials and Methods

##### Clinical Protocol

^{−1}· h

^{−1}(group 1, 2, 3, 4, or 5, respectively). The continuous infusions were administered by using a Graseby 3500 (Smiths Medical International, Watford, United Kingdom) infusion pump. Age (< 65 and ≥ 65 yr) and gender were stratified into all study groups.

*via*an 18- or 20-gauge catheter in a large forearm vein without any other fluids given simultaneously. As such, the remaining dead space was limited towards the dead space of the infusion catheter (less than 0.023 ml) filled with blood before starting the propofol administration. The propofol syringe and extension line were purged and pressurized before being connected to the infusion catheter. The moment the START button was pushed was considered as time 0.

*via*a facemask during the study period. When airway obstruction occurred, airway was secured gently without any stimulation.

##### Sample Acquisition, Handling, and Processing

*via*a closed-system sampling port at the arterial line (dead space 1.5 ml).

##### Measuring Propofol Cerebral Drug Effect

^{11}The BIS values were recorded every 5 s by Bispectrum Analyzer software developed by one of the authors

^{12}using an RS-232C connection between the BIS monitor and a laptop computer. The end of the study period was defined as 60 s after the end of last arterial sampling time.

##### Drug Assay

^{13}The lower limits of detection and quantification were 1.7 and 5.6 ng/ml, respectively. Measured Cp below the lower detection limit was not applied to the model development and measured Cp between the lower limit of detection and that of quantification was considered as the half value of the lower quantification limit.

##### Population Pharmacokinetic Modeling

##### Basic Model Structure.

_{tr}is a first-order transit rate constant describing the movement of propofol between TRANSIT compartments or TRANSIT compartment and central compartment in the TRANSIT model. The optimal number of TRANSIT compartments was assessed by stepwise addition or deletion of one TRANSIT compartment. All combinations of these models were compared.

_{1}, V

_{2}, V

_{3}, CL

_{1}, CL

_{2}, CL

_{3}, k

_{tr}, and LAG time; V being the distribution volume for compartment, CL being the clearance for the compartment) was estimated by using a log-normal distribution:

_{i}is the parameter value in the i-th patient, θ

_{TV}is the typical value of the parameter in the population, and η

_{i}is a random variable in the i-th patient with a mean of 0 and a variance of ω.

^{2}Interindividual variability is reported as ω, the SD of η in the log domain, which is approximately the coefficient of variation in the standard domain. Residual intraindividual variability was modeled by using constant coefficient of variance model and additive error model. The best model had the minimum values of Bayesian Information Criterion among the all evaluated model. These criteria were calculated by the following equations:

^{14}

##### Covariate Assessment for Pharmacokinetic Model.

*P*< 0.005, chi-square test with 1 degree of freedom) was considered significant. When the assessed model included an additional potential covariate without an increment of the number of the parameter estimates, a decrease in OFV was considered to obtain the potential covariate into the final model. Weight, height, age, and dose rate were assessed as potential covariates when the linear regression analysis indicated significant correlation (

*P*< 0.05) between a covariate and a parameter estimate. After significant regression, all potential covariates were tested on the parameters from the basic structural model (V

_{1}, V

_{2}, CL

_{1}, CL

_{2}, LAG, and k

_{tr}). Correlations between mean transit time and model parameters were also calculated to assess the accuracy of the inclusion of LAG or k

_{tr}into the model:

_{tr}is the number of transit compartments. For assessing the effects of a covariate on the model accuracy, we used not only OFV, but we also depicted the measured/predicted Cp and measured Cp

*versus*time. When this figure indicated severe bias in the assessed model, we excluded the model for further consideration.

##### Population Pharmacodynamic Modeling

##### Basic Model Structure.

_{e0}. A classic sigmoidal maximum possible drug effect (E

_{max}) model was used for describing the relation between propofol effect-site concentration (Ce) and the electroencephalographic measures of propofol with either measured plasma concentration (Cp) or predicted Cp. The time course of the predicted Cp was calculated during the pharmacodynamic modeling by using the individual patients’

*post hoc*Bayesian pharmacokinetic parameters estimated in the pharmacokinetic modeling. The Ce with measured Cp was calculated by using a “connect-the-dots” approach previously used by Schnider

*et al.*

^{15}and mathematically described previously.

^{6}The classic sigmoidal E

_{max}model is described by the following equation:

*e.g.*, the measured BIS value), E

_{0}is the baseline measurement when no drug is present, E

_{max}is the maximum possible drug effect, Ce is the calculated effect-site concentration of propofol, Ce

_{50}is the Ce associated with 50% maximum drug effect, and γ is the steepness of the concentration

*versus*response relation. The model parameters were estimated by using NONMEM. For Ce

_{50}, k

_{e0}, and γ, interindividual variability was permitted by using a log normal distribution (θ

_{i}= θ

_{TV}· e

^{- ηi}, which was the same equation in the pharmacokinetic model). Interindividual variability for E

_{0}was expressed by using a normal distribution: θ

_{i}= θ

_{TV}+ η

_{i}. Interindividual variability is reported as ω, the SD of η in the standard domain for E

_{0}. Residual intraindividual variability was modeled by using an additive error model.

_{max}model with a fixed E

_{max}at zero and an E

_{max}model with an estimated E

_{max}. Although one should be cautious with extrapolations, these two pharmacodynamic models were studied as one might criticize that the BIS might depict a maximum drug effect at a BIS value of zero. In addition, we included the BIS delay, defined as a time shift of applying the effect-site concentration to the sigmoidal E

_{max}model when BIS value and the effect-site concentration were applied to the model, which reflects the time delay for the BIS measurement. This delay is the result of the combination of averaging algorithms (

*i.e.*, smoothening rate), the delay caused by the BIS signal quality and the delay in adaptation of one of artifact rejection preprocessing steps.

^{6,16}In this study, the BIS delay was assessed as one value with no variation for all groups.

_{e0}to differ among groups.

##### Different K_{e0} Values for the Basic Model.

_{e0}for each basic structural model. When one k

_{e0}was determined for two or more groups, each k

_{e0}(different k

_{e0}) was used for each consecutive group (

*e.g.*, one k

_{e0}is for groups 1 and 2, and the other k

_{e0}is for groups 3, 4, and 5). The addition of k

_{e0}was considered statistically significant when the NONMEM objective function value (OFV or –2 log likelihood) decreased by at least 6.63 for an additional k

_{e0}(

*P*< 0.01, chi-square test with one degree of freedom).

##### Covariate Assessment for the Pharmacodynamic Model.

_{e0}, the final pharmacodynamic model was implemented to assess for possible covariates. In this process, a significant correlation (

*P*< 0.05) between a covariate and a parameter on the linear regression and a decrease in OFV by 6.63 (

*P*< 0.01, chi-square test with 1 degree of freedom) was an additional criterion for a covariate. Schnider

*et al.*reported that elderly patients are more sensitive to the hypnotic and electroencephalogram effects of propofol than younger persons.

^{15}Therefore, age was explored as a covariate for Ce

_{50}.

##### Model Evaluation.

*versus*time that shows the bias of the prediction

*versus*time.

*post hoc*Bayesian parameters) of the final model. For calculating these values, the following equations were used.

_{i}is the number of |PE| values obtained for the i-th subject.

#### Results

##### Pharmacokinetic Model Development

Table 2 Image Tools |
Fig. 3 Image Tools |
Table 3 Image Tools |

_{tr}. After the visual inspection of the correlation between dose rate and mean transit time (fig. 3), we decided to test these covariates for the final model. All covariates were examined for their potential linear relationship with the structural pharmacokinetic parameters (see table, Supplemental Digital Content 1, which lists OFVs for the examined models including potential covariates, http://links.lww.com/ALN/A548). The best model includes dose rate for k

_{tr}, age for LAG time, and weight for V

_{1}as covariates (fig. 1; see also table, Supplemental Digital Content 2, which shows the NONMEM representation of the final pharmacokinetic model, http://links.lww.com/ALN/A549). The parameter estimates in our final pharmacokinetic model are represented in table 3. Time courses of

*post hoc*Bayesian predicted concentration for all patients are shown on the right column in figure 2.

##### Model Evaluation for the Final Pharmacokinetic Model

*versus*both population and

*post hoc*Bayesian predicted concentrations (fig. 4A and 4B for the all groups; see also figures, Supplemental Digital Content 3, which include the goodness-of-fit plots for each group, http://links.lww.com/ALN/A550), and the time course of the measured/predicted Cp

*versus*time (fig. 4C–F for all groups; see also figures, Supplemental Digital Content 3, which include the time course for each group, http://links.lww.com/ALN/A550). These figures indicate that the model provided an accurate description of the observed data. The sampling point where measured and/or predicted Cp was below quantification limit of the drug assay was excluded from figure 4C–F (see figures, Supplemental Digital Content 4, which depict the detailed information for the below limit of quantification values, http://links.lww.com/ALN/A551). On the final population model, the predicted plasma concentrations showed larger variability in the first minute than that in the later time (fig. 4E).

*post hoc*prediction, respectively. These metrics indicated an acceptable bias and accuracy of the model.

##### Simulations with the Final Pharmacokinetic Model

_{tr}) also produces lag time (fig. 5B). For a typical patient, the influence of age and weight on the time course of the predicted arterial propofol concentration is plotted in figures 5D and 5F, respectively.

##### Pharmacodynamic Model Development

*post hoc*Bayesian individual predicted Cp for the pharmacodynamic examinations because NONMEM failed to calculate the converged estimates with measured Cp. Applying multiple k

_{e0}values into the model for different administration rates did not improve the model (see table, Supplemental Digital Content 5, which represents OFVs for the different k

_{e0}, http://links.lww.com/ALN/A552). Age was not included into the pharmacodynamic model because R square values between age and the

*post hoc*k

_{e0}value for the models including estimated E

_{max}or fixed E

_{max}at zero indicated no significance (R

^{2}= 0.001 or 0.039, respectively).

##### Final Pharmacodynamic Model

Table 4 Image Tools |
Table 5 Image Tools |
Fig. 6 Image Tools |

_{max}are shown in tables 4 and 5, respectively (see table, Supplemental Digital Content 6, which is the NONMEM representation of the final pharmacodynamic model, http://links.lww.com/ALN/A553). Interindividual variability is expressed as SD for E

_{0}and coefficient of variation for E

_{max}, Ce

_{50}, and γ. The k

_{e0}values were 0.414 min

^{−1}or 0.404 min

^{−1}for the model including estimated or fixed E

_{max}, respectively. The LAG times for BIS delay corresponding to the effect-site concentration were 19.7 or 20.3 s for the models including estimated or fixed E

_{max}, respectively. Figure 6 shows the measured BIS and

*post hoc*individual Ce

*versus*time, and the relation between the BIS values and the

*post hoc*individual predicted Ce of propofol.

_{max}model for the two final models indicates the similar relation between Ce and BIS when BIS value is more than 45 (fig. 7).

##### Pharmacodynamic Model Evaluation

*versus*time. Figure 8E and 8F depict the sigmoid curve and all measured BIS value

*versus post hoc*individual Ce. These figures indicate that the model provided an accurate description of the observed data.

_{max}or fixed E

_{max}at zero, respectively. Confidence intervals of 95% for parameter estimates are shown in tables 4 and 5 for the two final pharmacodynamic models. These results revealed that all the final estimated model parameters have acceptable typical values. The MADPE and MDPE values of the two final pharmacodynamic models are shown in table 6 for both population and

*post hoc*predictions. These metrics indicated an acceptable bias and accuracy of the model.

##### Correlation between Dose Rate and k_{e0}

*post hoc*k

_{e0}values

*versus*the dose rate of propofol in the final pharmacodynamic model. There were no significant differences among the groups using a Tukey multiple comparison test.

#### Discussion

##### Pharmacokinetics

^{6}The failure of these models is an expected consequence of a flawed assumption with conventional mammillary compartmental models of instantaneous mixing in the central compartment and immediate appearance of drug at the arterial sample side. More appropriate models are required to study the early phase kinetics, certainly in relation to possible covariates influencing these kinetics.

^{17}and lungs,

^{18}it was only our intention to model the delay by using a compartmental, mathematical approach. The LAG time and TRANSIT models were previously used to describe the absorption delay in orally given drugs.

^{8}Upton

*et al.*also used the concept of LAG time in their physiologically based model to describe the delay between the infusion site and the vascular mixing volume submodel and a TRANSIT model with three compartments in the lung submodel.

^{18}The recirculatory model by Avram and Krejcie also includes the chain of TRANSIT compartments for the central circulation in the model.

^{7}We assessed four basic structural models combining a conventional one-, two-, and three-compartmental model enlarged with the LAG time and/or TRANSIT model in an attempt to increase the accuracy of the traditional compartmental model for describing the time course of the propofol arterial concentration, especially during the first minutes of drug administration. We revealed that both LAG and TRANSIT compartments were required to describe our data accurately. We did not apply a physiologic model to depict the early drug disposition of propofol because this is much more complex and was not required to solve our hypothesis. However, we tried to add a recirculatory step to our compartmental model, similar to the work of Avram and Krecjie

^{7}; however, this did not result in a more accurate model when using our database (see figure, Supplemental Digital Content 7, which represents the examined compartmental model with recirculation-like characteristics, http://links.lww.com/ALN/A554).

^{19–21}our best basic model includes only two disposition compartments. This is because we focused on the first ten minutes to describe the early phase kinetics, especially in the first minutes. This was a study in patients, so it was considered unethical to draw more blood samples at a later time point. Thus, the time course of measured concentration of propofol excluded the information for γ-phase, where half-life is over 4 h.

^{22}We have to accept that we only studied the early disposition of propofol; as such, our pharmacokinetic model should not be extrapolated to predict the pharmacokinetic behavior of propofol at a later stage.

_{tr}was influenced by dose rate as shown in table 3 and depicted in figure 5. Several other covariates were tested and resulted in a final model, including age as a covariate for the LAG time and weight for V

_{1}(fig. 5). Other tested covariates didn’t improve our model.

*et al.*,

^{23}who reported that the time before propofol or indocyanine green emerges from the central vein was approximately 13 s, despite the difference of the sampling site and when taking into account our more sensitive lower limits of detections of propofol (50 or 1.7 ng/ml, respectively). The age distribution of our study population allowed us to test the influence of age on the propofol early phase kinetics. We found that the LAG time increased in older patients, as shown in figure 5. Although one may consider that lower cardiac output increased LAG time in older patients, there was no correlation between LAG time and cardiac output by using our data (see figures, Supplemental Digital Content 8, which show the relation between measured cardiac output and two model parameters of LAG time and transit rate constant, http://links.lww.com/ALN/A555). Others also found an influence of age on the propofol drug disposition, even when using classic multi-compartmental model without coverage of the early phase kinetics.

^{20}

_{1}of 0.0104 l/kg is smaller than observed by others.

^{19–21}When using conventional compartmental models, the central volume of distribution is dependent on the time of blood sampling and will expand with time. As such, with earlier blood sampling, one will obtain smaller central volumes. However, conventional compartmental models usually overestimate the central compartmental volume because they ignore the complexity of delay caused by intravascular mixing.

^{7}Weight was considered as a covariate for V

_{1}. As the covariate analysis is fully dependent on the study population, one should be careful when comparing models from propofol because different populations might result in different model-dependent covariates.

*versus*time revealed the initial large variability in the population. The model bias is within acceptable limits. In addition, very large variability is observed over the first minute after the start of propofol administration, which may be partially caused by the variability of transit time from infusion site to sampling site (fig. 4C and 4E). This result indicated that our final population model was not able to address the variability in the first minute, although the high-resolution arterial samples were taken. The individual

*post hoc*Bayesian propofol Cp are accurately described by the model. Figure 4 also depicts a good correlation between the individual

*post hoc*Bayesian predicted arterial concentration and measured arterial concentration. As shown in table 3, the 95% confidence intervals of the parameter estimates for our model calculated by the bootstrap analysis exclude zero, revealing the accuracy of the model. In addition, all parameter estimates calculated by NONMEM are similar to the median of the estimates calculated by bootstrap analysis (data are not shown). The MDAPE and MDPE demonstrate that the bias and the predictive inaccuracy are low for both population and individual

*post hoc*Bayesian estimates. A model is expected to be accurate if the typical value of the predictive inaccuracy is less than 25%.

^{21}Our results indicate that our final pharmacokinetic model has an acceptable accuracy and is in the range of previous publications.

^{21,24,25}

##### Pharmacodynamics

^{6,26}Our pharmacodynamics could be described by using a classic sigmoid E

_{max}model. As shown in figure 6, our propofol dosing strategy resulted in a significant decrease in BIS values in all groups. However, we didn’t obtain BIS levels lower than 25. As such, the typical value for E

_{max}in our model was 38.7. As one might criticize that we know that the BIS can go down to zero, we developed a second model fixing the E

_{max}at zero. Although we accept the danger of extrapolation, figure 7 shows that both models behaved similarly above a BIS value of 45. For both pharmacodynamic models, acceptable population and

*post hoc*Bayesian estimation were reached as proven by the accurate measured/predicted BIS values (fig. 8), bootstrap analysis, and MDAPE and MDPE (table 6).

^{6}To take care of these delays, we included a BIS delay into our modeling work. Under accurate signal quality conditions, the average delay can in theory be estimated around 10 and 15 s. We found a typical BIS delay of 19.7 s or 20.3 s when applying the model with the estimated or the fixed E

_{max}, which comes slightly longer than the theoretical delay. K

_{e0}values were also similar between two pharmacodynamic models (table 4 and 5). This similarity among the models illustrated a potential of extrapolation on the model including fixed E

_{max}at zero. However, one should be cautious with the possible lack of uniformity across the BIS between 100 and 0, as shown by Kreuer

*et al.*

^{27}

*et al.*approach for the pharmacodynamic modeling.

^{9}The two-step approach is preferred in pharmacokinetics-pharmacodynamics modeling because a large number of

*post hoc*parameter estimates causes the uncertainly to calculate the parameter estimates when one develops a pharmacokinetic and a pharmacodynamic model simultaneously.

_{e0}.

^{9,28}The standard model of k

_{e0}assumes that the rate of equilibration between the plasma and the site of drug effect is independent of the rate of drug administration; however, to the best of our knowledge, scientific validation for propofol is still incomplete. Doufas

*et al.*found that induction speed is not a determinant of propofol pharmacodynamics. In their study, they applied a combined pharmacokinetic-dynamic model, whereby the pharmacokinetics were described by using a classic pharmacokinetic model. Classic pharmacokinetic models do not describe the early phase kinetics accurately; therefore, the misspecification in the pharmacokinetic model over the first few minutes might have influenced their findings.

^{6}

_{e0}exists. As such, we aimed at studying the influence of the administration rate on the propofol plasma-effect site equilibration using the accurate individual

*post hoc*Bayesian predicted plasma concentration based on the enlarged pharmacokinetic model describing the early phase kinetics. We hereby focused on that range of infusion rates up to 1,200 ml/h available in clinically applied syringe pumps. Our final model applied one k

_{e0}value for administration rates between 10 and 160 mg · kg

^{− 1}· h

^{− 1}, hereby covering the infusion rates available in commercial syringe pumps. As such, our results suggest that a model with only one k

_{e0}is appropriate and can be applied during target-controlled infusion. Our conclusion that k

_{e0}is independent from dose rate is also confirmed in figure 9, which depicts the relation between dose rates and the patient individual k

_{e0}s. Of course, one must be careful to not extrapolate our conclusions to models including very fast bolus administration. In these manually administered bolus injections, multiple k

_{e0}s might be appropriate as found previously.

^{6}

_{max}model with BIS delay to describe the dynamics is accurate and takes care of the early phase pharmacology of propofol during infusion rate between 10 and 160 mg · kg

^{−1}· h

^{−1}. Kinetics and dynamics could be connected with a singular k

_{e0}not influenced by the propofol infusion rate. However, the early phase kinetics of propofol were influenced by the infusion rate; as such, dose rate was included as a covariate for k

_{tr}. LAG time was influenced by age. In addition, although our modeling was able to handle the early phase pharmacokinetics, the large variability of population predicted plasma concentration in the first 1 min after the start of propofol infusion suggested that there may be a limit to how much pharmacokinetic models can predict the plasma concentrations against large variability of plasma concentration.