Conventional compartmental pharmacokinetic models wrongly assume instantaneous drug mixing in the central compartment, resulting in a flawed prediction of drug disposition for the first minutes, and the flaw affects pharmacodynamic modeling. This study examined the influence of the administration rate and other covariates on early phase kinetics and dynamics of propofol by using the enlarged structural pharmacokinetic model.

Fifty patients were randomly assigned to one of five groups to receive 1.2 mg/kg propofol given with the rate of 10 to 160 mg . kg(-1). h(-1). Arterial blood samples were taken frequently, especially during the first minute. The authors compared four basic pharmacokinetic models by using presystemic compartments and the time shift of dosing, LAG time. They also examined a sigmoidal maximum possible drug effect pharmacodynamic model. Patient characteristics and dose rate were obtained to test the model structure.

Our final pharmacokinetic model includes two conventional compartments enlarged with a LAG time and six presystemic compartments and includes following covariates: dose rate for transit rate constant, age for LAG time, and weight for central distribution volume. However, the equilibration rate constant between central and effect compartments was not influenced by infusion rate.

This study found that a combined pharmacokinetic-dynamic model consisting of a two-compartmental model with a LAG time and presystemic compartments and a sigmoidal maximum possible drug effect model accurately described the early phase pharmacology of propofol during infusion rate between 10 and 160 mg . kg(-1). h(-1). The infusion rate has an influence on kinetics, but not dynamics. Age was a covariate for LAG time.

CONVENTIONAL multicompartmental mammillary pharmacokinetic models assume that drug added to the central compartment is instantaneously completely mixed and that this mixed plasma instantaneously appears in the arterial circulation. The failure of these models to accurately predict the time course of arterial drug concentration in the first minutes after administration is the result of this flawed assumption and has been extensively described.1–6Although 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,7More 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.8Others 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

A first-order process, described by k^{e0}, characterizes the drug transfer between the plasma and effect site.9When 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

Induction dose, time course of drug concentration, and drug effect of propofol are dependent on administration rate. Therefore, the inclusion of dose rate into kinetic-dynamic model in the early phase might optimize the description of the time course of early drug disposition and drug effect.6,10

The aim of this study was to investigate if a compartmental model including LAG and TRANSIT models can describe the time course of arterial propofol concentration in the first minutes after administration more accurately than a classic compartmental model. We also studied the influence of the pharmacokinetic model on the description of propofol’s cerebral drug effect as measured by the Bispectral Index (BIS). In addition, we aimed at investigating the influence of propofol’s administration rate on these early phase kinetics and dynamics.

## Materials and Methods

### Clinical Protocol

After securing local ethics committee (National Defense Medical College, Saitama, Japan) approval and written informed consent, we enrolled 50 patients who were American Society of Anesthesiologists physical status I or II, aged 30–76 yr, and scheduled to undergo elective gynecologic, orthopedic, or digestive surgery. Exclusion criteria included body mass index greater than 30, neurologic disorder, recent use of psychoactive medicine, and significant heart, hepatic, or renal impairment.

All patients were randomly assigned to one of five groups to receive 1.2 mg/kg propofol (1% Diprivan®; AstraZeneca K.K., Osaka, Japan) given as continuous infusion with the rate of 160, 80, 40, 20, or 10 mg · kg^{−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.

Propofol was infused directly *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.

No other drugs were given during the study period. All patients received 100% oxygen *via* a facemask during the study period. When airway obstruction occurred, airway was secured gently without any stimulation.

### Sample Acquisition, Handling, and Processing

A 22-gauge arterial catheter was inserted in the radial artery at the opposite site of the venous catheter for propofol infusion.

To ensure accurate blood sample quality, a continuous blood flow during the study period of around 50 ml/min was generated by using a negative pressure blood drawing device (Hemoquick®; Terumo, Tokyo, japan). Arterial blood samples (1 ml each) were drawn by using heparine-coated 2.5-ml syringes *via* a closed-system sampling port at the arterial line (dead space 1.5 ml).

Arterial blood samples were taken at 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, and 60 s during the first 1 min in all groups. Samples were then collected every 10 s until 120, 150, or 200 s in groups 1, 2, and 3, respectively. In group 4, samples were taken at 70, 80, 90, 120, 150, 180, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300, and 310 s, and in group 5 at 70, 80, 90, 120, 150, 180, 210, 240, 300, 360, 420, 430, 440, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, and 530 s. These blood samples were centrifuged within an hour of collection. The plasma was transferred to polyethylene tube and kept at –20°C until assayed.

### Measuring Propofol Cerebral Drug Effect

Propofol cerebral drug effect was monitored by using the BIS, which was derived from the frontal electroencephalogram and calculated by the A-2000 BIS monitor (version 4.0; Aspect Medical Systems, Norwood, MA) by using the three BIS Sensor electrodes. The smoothening rate of the BIS monitor was set at 15 s. BIS values range from 100 to 0, with lower values denoting more drug effect.11The BIS values were recorded every 5 s by Bispectrum Analyzer software developed by one of the authors12using 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

Arterial plasma concentrations (Cp) of propofol were determined by high-performance liquid chromatography (RF-550, CTO-10AS, LC-10AD, SIL-10AD, SCL-10A, and DGU-14A; Shimadsu, Kyoto, Japan) with fluorescence detection at 310 nm after excitation at 276 nm.13The 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.

The following structural models were evaluated; conventional 1-, 2-, and 3-compartmental models, conventional one-, two-, and three-compartmental models enlarged with LAG model, conventional one-, two-, and three-compartmental models enlarged with a TRANSIT models, and conventional one-, two-, and three-compartmental models enlarged with a LAG and TRANSIT model (fig. 1). All compartmental models had a central compartment, and one or two peripheral compartments, if necessary. The LAG time and/or the TRANSIT model describes the time required for the venously injected drug to reach the sampling site of arterial blood. The k^{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.

All model parameters were estimated by using NONMEM VI (GloboMax LLC, Hanover, MD). The interindividual variability on each basic model parameter (V^{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:

where θ^{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 ω.2Interindividual 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:

where OFV is –2 log likelihood calculated as the objective function value (OFV) by NONMEM, K is the number of parameter estimates in the model, and N is the sample size of the data set. The Bayesian Information Criterion value was used to select the best structural model because the Bayesian Information Criterion would provide a more appropriate selection than the OFV when the different structural models are compared.14

#### Covariate Assessment for Pharmacokinetic Model.

The best basic structural model was implemented to assess the covariates by using a stepwise forward addition and backward elimination approach. In this process, a decrease in OFV by 7.88 (*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:

where N^{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.

The effect-site was assumed to be linked to the plasma by a compartment with a first-order equilibrium constant of k^{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.* 15and mathematically described previously.6The classic sigmoidal E^{max}model is described by the following equation:

where Effect is the electroencephalographic effect (*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.

In the pharmacodynamic modeling, we compared an E^{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,16In this study, the BIS delay was assessed as one value with no variation for all groups.

For the pharmacodynamic model approach, we assumed that the underlying sigmoidal model that describes the relation between Ce and the drug effect is not influenced by the dose rate of drug administration. This is a fundamental assumption underlying the standard model of the effect-site. Therefore, we concurrently estimated the model parameters for all five groups, only permitting the value of k^{e0}to differ among groups.

#### Different K^{e0}Values for the Basic Model.

We assessed the different k^{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.

After the assessment of the different k^{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.15Therefore, age was explored as a covariate for Ce^{50}.

#### Model Evaluation.

For a basic internal evaluation of the final pharmacokinetic model, a goodness-of-fit plot was used. The plot included the correlation between individual predicted concentrations by the Bayesian estimates and measured concentrations of propofol. We also depicted measured/predicted Cp or measured/predicted BIS *versus* time that shows the bias of the prediction *versus* time.

Bootstrap analysis was used for advanced internal evaluation of both pharmacokinetic and dynamic final model. One thousand bootstrap resampling data sets for both were created by sampling data from the original data set with replacement. Sample size of resampling data sets was the same as that of original data set. The final pharmacokinetic and dynamic model, including all significant covariates, was fitted to each of the resulting data sets. Confidence intervals of 95% were obtained for each parameter estimate as the 25th and 97.5th percentiles.

Median absolute prediction error (MDAPE) and median prediction error (MDPE) were also calculated as the metrics of the accuracy by using the individual parameters (population or *post hoc* Bayesian parameters) of the final model. For calculating these values, the following equations were used.

Percentage prediction error (PE) of the predicted CP:

Prediction error (PE) of the predicted BIS:

PE is an indication of the bias of the achieved concentrations, and the absolute value PE (|PE|) is a measure of the precision (inaccuracy).

MDAPE indicates the inaccuracy of the Cp or BIS predictions in the i-th subject:

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

MDPE reflects the bias of the Cp or BIS predictions in the i-th subject:

To illustrate the final pharmacokinetic model, especially for the effect of covariates, several simulations were implemented.

## Results

### Pharmacokinetic Model Development

All recorded data were used. The airway was secured with gentle jaw thrust in two patients, one in the group 4 and the other in group 5, because of airway obstruction. No patient experienced hemodynamic instability during the study. The demographics for the patients in the five groups are shown in table 1.

The left column in figure 2shows the time course of the measured Cp of propofol for the five groups, respectively. There were no samples taken after 280 s for one patient in group 4 and after 460 s for another patient in group 5 due to a protocol violation. One sample at 70 s in another patient in group 4 was not available because of a technical problem with the sample line. No samples revealed concentrations between the lower limit of detection and that of quantification.

For the basic structural model, a two-compartmental model enlarged with a LAG and TRANSIT model including six TRANSIT compartments (fig. 1) was selected as the best basic model (table 2). Linear regression analysis indicated several potential covariates, including dose rate for LAG and k^{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

Figure 4represents the goodness-of-fit plot for the measured plasma concentrations *versus* both population and *post hoc* Bayesian predicted concentrations (fig. 4A and 4Bfor 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–Ffor 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).

One thousand bootstrap runs were obtained with the same sample point at the calculating estimates of the final model. The number of successfully converged resamples was 798. Confidence intervals of 95% for parameter estimates are shown in table 3. The intervals for thetas exclude zero values. These results revealed that all the final estimated model parameters have acceptable typical values. MDAPE and MDPE of the final pharmacokinetic model were 23.6% and –1.4%, and 8.9% and 0.8% for population and *post hoc* prediction, respectively. These metrics indicated an acceptable bias and accuracy of the model.

### Simulations with the Final Pharmacokinetic Model

We performed several simulations to show the influence of the covariates applied our final model as shown in figure 5. The influence of dose rate on the initial time course of the propofol arterial concentration as predicted by the final model is shown in figure 5Bfor a typical patient. Lower infusion rate results in longer time between the start of infusion and the start of increasing predicted concentration than faster infusion rate as transit rate constant (k^{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

Five patients (one each in groups 1, 2, and 4 and two in group 5) were excluded because of technical problems with the BIS recordings (automatic recheck of the BIS monitor during the study period due to insufficient sensor fitting). Therefore, 45 patients were included into the pharmacodynamic analysis.

After initial attempts, we decided to use the *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

The parameter estimates in our final pharmacodynamic model including estimated or fixed E^{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 6shows 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.

The sigmoid curves of the sigmoidal E^{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

For both pharmacodynamic models, figures 8A–8Dshow the goodness-of-fit plot for the time course of the measured/predicted BIS *versus* time. Figure 8E and 8Fdepict 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.

One thousand bootstrap runs were obtained with the same sample point at the calculating estimates of the two final models. The numbers of successfully converged resamples were 891 or 904 for the models including estimated E^{max}or fixed E^{max}at zero, respectively. Confidence intervals of 95% for parameter estimates are shown in tables 4 and 5for 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 6for both population and *post hoc* predictions. These metrics indicated an acceptable bias and accuracy of the model.

### Correlation between Dose Rate and k^{e0}

Figure 9shows the *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

We found that a multi-compartmental model including a LAG time and TRANSIT model described the time course of arterial propofol concentration in the first minutes after administration more accurately than a classic compartmental model. In addition, we found that the propofol administration rate influenced the early phase kinetics but not the dynamics. Age and weight were also defined as covariates in our pharmacokinetic model.

### Pharmacokinetics

Conventional mammillary multi-compartment pharmacokinetic models for propofol are well accepted and form the basis for clinically applied target-controlled infusion systems. Unfortunately, these classic models insufficiently predict the time course of the arterial propofol drug concentration for the first minutes of drug administration.6The 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.

Before being able to analyze the influence of covariates such as the infusion rate on the dose-response of propofol, we needed the model including the delay between drug administration and appearance of the drug at the site of arterial blood sampling. Although we realized that this delay might be explained mechanistically or by physiologic phenomena such as transit delay due to passage through the forearm vein17and lungs,18it 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.8Upton *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.18The recirculatory model by Avram and Krejcie also includes the chain of TRANSIT compartments for the central circulation in the model.7We 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 Krecjie7; 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).

The final structural model was a two-compartmental model enlarged with a LAG time and TRANSIT model. Although most previously published conventional compartmental pharmacokinetic models of propofol have three compartments,19–21our 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.22We 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.

Covariate analysis started by investigating the influence of dose rate on the early phase kinetics. K^{tr}was influenced by dose rate as shown in table 3and 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.

These findings are consistent with the results from He *et al.* ,23who 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

As shown in table 3, V^{1}of 0.0104 l/kg is smaller than observed by others.19–21When 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.7Weight 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.

We tested the internal validity of our model by using several analyses. As shown in figure 4, the population measured/predicted Cp *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 4also 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%.21Our results indicate that our final pharmacokinetic model has an acceptable accuracy and is in the range of previous publications.21,24,25

### Pharmacodynamics

The propofol cerebral drug effect was measured by using BIS, as applied in various previous studies.6,26Our 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 7shows 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).

When analyzing the raw data, a certain delay in the BIS behavior might be observed. This delay may reflect the combination of the averaging algorithm to calculate the BIS and the delay in adaptation of one of the artifact rejection processing steps.6To 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

To develop the pharmacodynamic model, we used the traditional two-step described by Sheiner *et al.* approach for the pharmacodynamic modeling.9The 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.

Propofol transfer between the plasma and effect-site can be modeled as a first-order process characterized by k^{e0}.9,28The 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

We studied whether a real influence of the propofol administration rate on k^{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

In conclusion, we found that a combined pharmacokinetic-dynamic model consisted of a two-compartmental model with a LAG time and TRANSIT compartments to describe the kinetics and a sigmoid E^{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.

The authors thank Isao Fukuda, M.D., Ph.D. (Staff Anesthesiologist, Department of Anesthesiology, Tokyo Hospital, Tokyo, Japan) for providing the methodology for arterial blood sampling.

## References

*in vitro*haemolytic potential of a solution of a new medicine. Comp Haematol Int 1996; 6:35–41