Previous studies have reported conflicting findings concerning the influence of age on the pharmacokinetics of fentanyl, alfentanil, and sufentanil. Lemmens and associates 
considered that a possible factor explaining the differences between various studies could be the effect of gender. Conflicting pharmacodynamic findings for the opioids have also been reported. Scott and Stanski 
found that brain sensitivity to opioids, as determined by changes in the spectral edge (SE sub 95) of the electroencephalogram (EEG), increased significantly with age. In contrast to these findings, Lemmens and coworkers [3,4]
could not demonstrate an effect of age on the plasma concentration‐effect relations for perioperative stimuli.
The aim of this study was to determine the influence of subject covariates (particularly age and gender) on the pharmacokinetics and pharmacodynamics of the new short‐acting opioid remifentanil, and to assess the interindividual variability for the age range of 20 to 85 y. Using the spectral edge as a measure of drug effect, we investigated the pharmacokinetics and pharmacodynamics of remifentanil in 65 persons stratified by age and gender, and we prospectively applied the derived models in another 15 persons ages 41 to 84 y.
Materials and Methods
After obtaining Institutional Review Board approval and informed consent, 65 volunteers classified as American Society of Anesthesiology (ASA) physical status 1 and 2 were enrolled in the trial. Our goal in volunteer recruitment was to have nearly equal representation in three age groups: young (20–40 y), middle‐aged (40–65 y), and elderly (> 65 y). The ages ranged from 20–85 y. Thirty‐eight men and 27 women were studied. Remifentanil blood concentration and EEG effect data were collected in four separate phases. In the first phase, 10 healthy men classified as ASA physical status 1 whose ages ranged from 22 to 38 y received a constant‐rate infusion of remifentanil at doses ranging from 1 to 8 micro gram [centered dot] kg sup ‐1 [centered dot] min sup ‐1 for 20 min. 
In the next phase, ten more healthy men (ages 23 to 39 y) were enrolled in a cross‐over study to investigate the comparative pharmacokinetics and pharmacodynamics of alfentanil and remifentanil. 
In the remifentanil arm of the trial, the volunteers received an infusion of 2 or 3 micro gram [centered dot] kg sup ‐1 [centered dot] min sup ‐1 for 10 min. 
In a third phase, ten healthy women volunteers (ages 20 to 33 y) were studied according to the same study protocol. In a fourth phase, 35 healthy men and women (ages 40 to 85 y) were enrolled. All volunteers in this phase received a remifentanil infusion at 3 micro gram [centered dot] kg sup ‐1 [centered dot] min sup ‐1 until there was no further change in the EEG, as measured by the spectral edge frequency (infusion duration, 4–15 min). We studied 15 more persons using the same protocol (11 men and 4 women between the ages of 41 and 84 y) from whom we collected EEG but not concentration data.
Sample Acquisition, Handling, and Processing
Remifentanil was administered through an 18‐ or 16‐G intravenous catheter inserted into a forearm vein. A 20‐G arterial catheter was inserted into the radial artery for sample collection. During all phases, frequent 3‐ml arterial blood samples were taken at preset intervals until 4 h after the end of the infusion. The rapid metabolism of remifentanil by nonspecific blood and tissue esterases required that all esterases be inactivated immediately after the arterial samples were taken. This was achieved by adding acetonitrile to the sample to denature the plasma proteins. Subsequently, methylene chloride was added to the mixture to extract the remifentanil into the organic phase, which was then separated and stored at ‐70 degrees Celsius until assay. The remifentanil analyses were performed using a high‐resolution, gas chromatographic‐mass spectrometry assay with a quantitation limit of 0.1 ng/ml and a paired aliquot coefficient of variation of less than 15% for concentrations greater than 0.1 ng/ml. [7,8]
An initial exploratory two‐stage pharmacokinetic analysis was performed to determine whether a two‐ or three‐compartment mammillary model best described each individual's data. To develop the final population pharmacokinetic model, we used the approach of Mandema and colleagues. 
Two‐ and three‐compartment population pharmacokinetic models without covariates were estimated with NONMEM* using the “first order conditional estimation” method and “eta‐epsilon interaction” to reduce the influence of model misspecification. The interindividual error on each of the model parameters (V1
, Cl sub 2, and Cl3
) was modeled using a log‐normal variance model: Equation 1
is the value of the parameter in the individual, thetaTV
is the typical value of the parameter in the population, and eta is a random variable with mean zero and variance omega sup 2. Technically, omega is the standard deviation of eta in the log domain. However, when omega is small, it is approximately the coefficient of variation (%CV) of P. To simplify the interpretation of our results, we report the variability of the pharmacokinetic parameters, P, as “%CV,” understanding that this is an approximate interpretation of the modeled variability. A “constant c.v.” model was used for the residual intraindividual error.
Empirical Bayesian [10,11]
estimates of the individual pharmacokinetic parameters were obtained based on the typical values and variances of the estimated volumes and clearances. Unit disposition functions were calculated from these Bayesian estimates and plotted between the times of the first and last observations. To confirm that the population parameters provided a good description of the individual unit disposition functions, we superimposed the unit disposition function for the typical person on this graph. The median of the weighted residuals (calculated as the difference between the individual and the population unit disposition functions at the time of each sample) was computed as a measure of the bias over the first 10 min.
A generalized additive model (GAM) was used to permit identification of linear and nonlinear relations between the Bayesian estimates of individual volumes and clearances and the volunteer covariates. We provided the GAM function (as implemented in S‐PLUS**) with a list of covariates available for analysis. The covariates analyzed were age, gender, weight, height, body surface area, and lean body mass (LBM). Body surface area was calculated from the weight (measured in kilograms) and height (measured in centimeters) using the formula of DuBois and DuBois 
: Equation 2
. Lean body mass was calculated from the gender, weight (in kilograms), and height (in centimeters)***: Equation 3
, Equation 4
For each of the covariates in this list, we provided the GAM function with an ordered regimen of possible candidate forms describing how each covariate could enter the model. We allowed three forms for each covariate. Each covariate could (1) be omitted from the model; (2) enter the model linearly; or (3) enter the model nonlinearly. We chose a natural cubic spline with knots at the 33rd and 66th percentile of the data as a flexible function to search for potential nonlinear relations. A nonlinear form was not included for gender. The GAM function performed a stepwise search to find the significant covariates for each of the volumes and clearances and to find the best form (linear or nonlinear) of each important covariate. Not all possible combinations of covariates and forms were examined. The model that resulted in the biggest decrease in the Akaike Information Criterion was returned by the GAM function as the best model. Ninety‐five models were examined in this step. The GAM function is described in more detail by Hastie. 
The structural model incorporating covariates identified by the GAM analysis was then evaluated with NONMEM to develop the final parameter estimates. The final model was examined for parsimony by excluding individual covariates and demonstrating a statistically significant increase in the NONMEM objective function, and by analysis of the standard errors of the estimated parameters.
As previously defined, we used the weighted residual to describe the quality of the prediction in each individual, calculated as (measured‐predicted)/predicted. 
We examined the performance of the complex pharmacokinetic model in the worst case for young, middle‐aged, and elderly subjects. These persons were selected based on the worst mean absolute weighted residual in each age group. We also examined the quality of the prediction of the pharmacokinetic models for the population by calculating the median absolute weighted residual (MDAWR) and the mean of the individual mean absolute weighted residuals (MAWR) as described by Kataria and associates, 
and by examining the pattern of the weighted residuals, expressed as measured/predicted value. 
In anticipation of developing dosing guidelines based on weight, we examined the performance of a weight proportional model and a lean body mass proportional model, in which all of the volumes and clearances are assumed to increase in proportion to the individual's weight or lean body mass, respectively.
The EEG was digitized and the spectral edge calculated as previously reported. 
We performed an initial two‐stage pharmacodynamic analysis, calculating the effect site concentration using the individual two‐ or three‐compartment pharmacokinetic parameters obtained during the two‐stage pharmacokinetic analysis. To develop the final population pharmacodynamic model, we again used the approach of Mandema and colleagues. 
NONMEM's first‐order conditional estimation method was used to analyze the pharmacodynamic data (spectral edge frequency) with a sigmoid Emax
model. Equation 5
is the measured baseline effect, Emax
is the maximum measured effect, Ce is the apparent effect site concentration, EC50
is the concentration that causes 50% of the maximum effect, and gamma describes the steepness of the concentration response relation. The apparent effect site concentration was calculated for each individual based on the individual dosing information, previous Bayesian estimates of the individual pharmacokinetic parameters (which entered into the pharmacodynamic model as known parameters or “fixed effects”), and the equilibration rate constant, ke0
(estimated by NONMEM). We used a log‐normal variance model to describe the residual interindividual error on each of the model parameters (ke0
, gamma) and an additive model to describe the residual intraindividual error (in contrast to the constant c.v. model used for the pharmacokinetic analysis). The relation between the Bayesian estimates of the individual pharmacodynamic parameters and the volunteer covariates was also examined using a GAM, as described previously. Seventy models were evaluated in this step.
As described for the pharmacokinetic analysis, the pharmacodynamic model incorporating covariates identified by the GAM analysis was evaluated with NONMEM to develop the final pharmacodynamic model. The final model was examined for parsimony by excluding individual covariates and demonstrating a statistically significant increase in the NONMEM objective function, and by analyzing the standard errors of the estimated parameters.
In contrast to the kinetic analysis, we used the unweighted residual (observed‐predicted) to describe the quality of the prediction in each individual. We examined the performance of the simple and complex pharmacodynamic model in the worst case for young, middle‐aged, and elderly persons. We also examined the quality of the prediction of the pharmacodynamic models for the population by calculating the median absolute residual and the mean of the individual mean absolute residuals (MAR). To provide some sense of the magnitude of these values, we also reported them as a percentage of the pharmacodynamic range (E0 ‐Emax) of the respective model.
We used the dose‐versus‐response data from the 15 additional participants to compare prospectively the performance of the simple (without covariates) and complex (with covariates) pharmacokinetic and pharmacodynamic models. We examined the predictive performance of the simple and complex pharmacokinetic and pharmacodynamic models in the the worst case for the middle‐aged and elderly subjects. We examined the quality of the predictions of the combined model by calculating the median absolute prediction error and the mean of the individual mean absolute prediction errors (MAPE) as measures of the precision of the models. We also examined the median prediction error and the mean of the individual mean prediction errors (MPE) as measures of the bias of the models, as detailed by Varvel and coworkers. 
However, we defined prediction error as the difference between the observed EEG effect and the predicted EEG effect, in keeping with our use of an additive variance model for the residual intrasubject error. These values were reported as a percentage of the pharmacodynamic range (E0
) of the respective model.
All 65 participants enrolled completed the study. Table 1
summarizes their distribution by age group and gender, and Figure 1
shows the individual values by age, gender, weight, height, LBM, and body surface area.
The initial two‐stage pharmacokinetic analysis found that the concentration‐time profiles for approximately one half of the individuals were best described using a two‐compartment model, and the other half were best described using a three‐compartment model. The plots of V1, Cl sub 1, and Vdas (calculated as V1 + V2 for the two‐compartment models, and as V1 + V2 + V3 for the three‐compartment models) showed decreasing trends with increasing age.
The population model derived using NONMEM found that the pharmacokinetics of remifentanil were best described by a three‐compartment model (P << 0.01). The parameter estimates and performance of the three‐compartment model without covariates are shown in Table 2
, together with the results of the weight proportional and the LBM proportional models. The LBM proportional model resulted in a greater improvement in the objective function than the weight proportional model, although this was associated with only a small improvement in the MDAWR and (MAWR) with bar. The main improvement resulting from the LBM proportional model was a 10% decrease in the %CV of V1
, together with small decreases in the %CV of V2
) shows the unit disposition functions obtained with the simple three‐compartment model. The unit disposition functions were calculated from the time of the first to the time of the last detectable concentration. Although sampling was performed for 4 h beyond the end of infusion, the remifentanil concentrations in many individuals were less than the limit of detection before the third exponential phase was apparent. The simple model (without covariates) described the early clinically important component of the unit disposition functions well, with a small negative bias (‐5.75%) for the first 10 min.
The GAM method was used to identify potentially significant covariates. The best models included age for all volumes and clearances. In addition, LBM was a covariate of V1
, and Cl1
in the best models. These were the same parameters for which inclusion of LBM as a covariate decreased the %CV in the LBM‐adjusted model described before. No significant nonlinear relations were detected. Figure 3
shows the individual Bayesian estimates of V1
as a function of age.
An initial model with age on all parameters and LBM on V1
, and Cl1
was implemented in NONMEM. The covariates were modeled around their median value, such that the intercept term reflects the typical value for a person of median age and LBM. This initial model was further refined by deleting certain parameters according to the criteria outlined by Mandema and colleagues. 
Only the linear relations between age and V3
did not remain in the final model (Table 3
). From the age of 20 to 85 yr, we found that typical values for V1
decreased by approximately 25% and 33%, respectively.
The more complex model accounted for a 9% to 16% decrease in the %CV of the structural model parameters. In addition, the standard errors on the variance of eta for all parameters except V3
were decreased by 15% (V1
) to 50% (Cl1
). This more complex pharmacokinetic model showed a substantial improvement in ‐2 log likelihood (213.5; P < 0.01) and was associated with an improvement in the MDAWR and MAWR with bar of approximately 5% and 4%, respectively. The improvement was consistent across all age groups, as assessed by the worst mean absolute weighted residuals. Figure 4
shows the weighted residuals for the simple model, the LBM proportional model, and the complex pharmacokinetic model. The improvement with the addition of both LBM and age as significant covariates is seen in the lower panel.
Four participants were excluded from the pharmacodynamic analysis because of very noisy or no EEG data. The initial two‐stage pharmacodynamic analysis detected an age‐related decrease in EC50
. Figure 5
shows the individual Bayesian estimates of the concentration‐response relations, together with the typical response of the population model without covariates. This figure suggests that there is considerable interindividual variability in E0
, and gamma. Table 4
shows the parameters and performance of the simple pharmacodynamic model. The median absolute residual for the simple model was 2.11 Hz (15% of the dynamic range).
The GAM method was used to identify potentially significant covariates. The best models included an age‐related linear decrease in k sub e0 and EC50
. No significant nonlinear relations were detected. The decreases in ke0
with age were confirmed by subsequent analysis with NONMEM. Table 5
shows the final age‐adjusted pharmacokinetic model. The parameters were modeled about the median age, such that the intercept reflects the typical value for a 40 y old. Figure 6
illustrates the doubling of t1
and the halving of EC sub 50 that occur from ages 20 to 85 y in the typical individual. The complex pharmacodynamic model showed a significant improvement in ‐2 log likelihood (27; P < 0.01) associated with only a small improvement in the median absolute residual and MAR with bar. This model resulted in an 8% decrease in the %CV of the ke0
and a 9% decrease in the %CV of the EC sub 50. As with the pharmacokinetic analysis, the standard errors on the estimates of the variance of eta were decreased for those parameters with significant covariates in the complex model. The standard errors of ke0
were decreased by 17% and 11%, respectively. The worst median absolute variance in the young and middle‐aged groups also improved. Some improvement in the pattern of the residuals was noted after infusion (Figure 7
The 15 participants (11 men, 4 women) in whom the simple and complex pharmacokinetic and pharmacodynamic models were tested prospectively were ages 41 to 84 y and had LBMs of 41 to 74 kg. The median absolute prediction error of the simple model was 23% of the pharmacodynamic (E0 ‐Emax) range. The median absolute prediction error of the complex model was 19% of the pharmacodynamic range. The improvement in median absolute prediction error in this prospective analysis suggests that the inclusion of age and LBM in the complex model will allow more accurate model performance. However, this improvement came at the cost of an increase in bias: The median prediction error was 0.5% with the simple model and 7% with the complex model.
Unique features of remifentanil are its rapid clearance and rapid ke0
, resulting in a rapid onset and offset of drug effect. It is tempting to speculate that these characteristics will make remifentanil an easy drug to titrate, and that clinicians will not need to consider patient covariates when choosing a dosing regimen. However, the rapid onset of drug effect may be accompanied by rapid onset of adverse events such as apnea and muscle rigidity. The rapid offset of drug effect can cause severe patient pain when the anesthesiologist is ill equipped to deal the problem, such as when the patient is being transported to the recovery room. In a companion manuscript, 
we explore the dosing implications of the influence or age, lean body mass and intersubject variability on remifentanil dosing guidelines.
Previous studies have reported conflicting findings concerning the influence of age and gender on the pharmacokinetics of opioids. Helmers and associates 
compared alfentanil kinetics between young and elderly patients after a single bolus administration. They found that the plasma clearance of alfentanil is less rapid and the elimination half‐life is longer in elderly patients than in young ones. In contrast, Scott and Stanski 
found no age‐related changes in the pharmacokinetic parameters of fentanyl or alfentanil. Hudson and colleagues 
found that patient age was positively correlated with both alfentanil Vdas
and elimination half‐life. Lemmens and associates 
considered that the effect of gender could be a possible factor to explain the differences among various studies. These authors studied the effects of age on the pharmacokinetics of alfentanil in female and male patients undergoing lower abdominal surgery with nitrous oxide‐alfentanil anesthesia. They found a significant negative correlation between plasma alfentanil clearance and age in women but not in men and concluded that the effects of age on the pharmacokinetics of alfentanil depend on gender. Matteo and coworkers 
found that the initial volume of sufentanil distribution was significantly smaller in elderly patients, but that the elimination half‐lives, plasma clearances, and total volumes of distribution were similar for elderly and younger participants. These authors believed that age‐related differences in the action of sufentanil could not be explained by the observed differences in the initial volume of distribution and that alterations in pharmacodynamics appear to be more important in the prolonged opioid effect seen in the elderly. Lehmann and colleagues 
found no significant correlation among sufentanil's half‐life, the volume of distribution, or clearance with age or body weight. In a more recent study, Helmers and colleagues 
reported that the sufentanil plasma concentration‐time curves for young adult and elderly participants could be superimposed. The elimination half‐life, clearance, and volume of distribution did not differ between both groups. They concluded that the age‐related differences in the action of sufentanil cannot be explained by its pharmacokinetic properties.
Conflicting pharmacodynamic findings for the opioids have also been reported. Scott and Stanski 
found a 50% decrease in the dose requirement of fentanyl and alfentanil with increasing age. They found that brain sensitivity, as determined by EEG changes, increased significantly with age. Lemmens and coworkers [3,4]
investigated the effects of age on the pharmacodynamics of alfentanil, as a supplement to nitrous oxide anesthesia, in 14 younger and 14 older women undergoing curative surgery for primary breast cancer. They could not demonstrate an effect of age on the plasma concentration‐effect relations for response to intubation, skin incision, and the probability of needing naloxone at the end of surgery. Ausems and colleagues 
examined the response of 37 patients between 23 and 54 y to perioperative stimuli using alfentanil and 66% nitrous oxide in oxygen. For the single events of intubation, skin incision, skin closure, and postoperative ventilation, each patient's clinical state was categorized as “responsive” or “nonresponsive” based on predefined criteria. The data from all patients were pooled and analyzed with logistic regression to define alfentanil concentrations required along with 66% nitrous oxide to provide a satisfactory state of general anesthesia. They also recorded multiple observations during surgery to determine the quantal response to stimulation in each patient. No relation between pharmacodynamic parameters and age were reported.
We evaluated several different pharmacokinetic models for the following reasons. The two‐stage pharmacokinetic analysis identified both two‐ and three‐compartment models, which made it difficult to average the parameters to obtain a two‐stage population parameter set. In addition, one half of the individuals did not have the third‐compartment pharmacokinetic parameter estimates necessary for the covariate analysis. The mixed‐effects modeling approach resolved these difficulties. Bayesian estimation of individual pharmacokinetic and pharmacodynamic parameters resulted in shrinkage of the parameter estimates toward the population mean characteristic of population analysis while still providing the individual estimates necessary for the covariate analysis.
The weight‐proportional models may offer advantages to the clinician because they enable dosing guidelines to be derived on a per‐kilogram basis. However, both simple and weight‐adjusted models must be evaluated before arbitrarily selecting per‐kilogram dosing as the preferred dosing model. Lean body mass entered the final complex pharmacokinetic model in a time‐dependent and nonlinear manner, and its influence was less important than age. The complex models enable us to develop dosing guidelines for clinicians that will account for the estimates of interindividual variability in pharmacokinetic and pharmacodynamic parameters. The complex models are probably the best choice for computer‐controlled drug delivery, although they are also appropriate for developing dosing nomograms and specific guidelines for use in special populations.
Our study involved more participants than did the opioid studies we cited. Because many covariates were correlated (such as gender, height, weight, LBM, and body surface area), and the data set was large, an efficient search strategy was needed to relate the participant covariates to pharmacokinetic and pharmacodynamic parameters. The GAM procedure provided a method to efficiently identify significant covariates. It was also able to build a model in the presence of correlated covariates, which presented a problem for more graphical methods of covariate analysis. 
For instance, in each age group, nearly all the men had a greater LBM than did the women (Figure 1
). Simply plotting the individual estimates of V1
against gender revealed that V1
was significantly less in the women. Although gender enters into the calculation of LBM, the GAM procedure could detect that LBM was a more important covariate than gender. Even with this large data set, the GAM procedure was very fast. The entire GAM analysis for the pharmacokinetic and pharmacodynamic data was complete in approximately 2 min. In contrast, evaluation of the final models with NONMEM required several months of continuous processing on a 90‐MHz Intel Pentium computer processor running Windows NT (Microsoft Corp., Redmond, WA) and an optimized 32‐bit version of NONMEM (Intel, Santa Clara, CA). An additional benefit of the GAM procedure was that it provided good initial estimates for the NONMEM analysis. However, we do not claim that by using the GAM procedure and the subsequent testing of the models with NONMEM we have necessarily found the optimal model. There may be other models that we have not tested that describe the observations better than the models we described. Testing every model with NONMEM would require years of effort with the available resources, and the likely improvement (assuming better models exist) would be small given the good performance of the models reported. Therefore, we restricted our exploration to those models that were identified as promising by the GAM analysis.
In a smaller study of 24 patients, Westmoreland and coworkers 
concluded that age did not influence the total clearance of remifentanil. In contrast, we found that age is an important covariate in the pharmacokinetics of remifentanil. We also found a reduced central volume of distribution in the elderly, which will result in higher initial blood concentrations after a given bolus dose than in younger persons.
Although the canonical univariate parameter for opioids reported by Gregg and colleagues 
and Gambus and associates 
is a more consistent measure of opioid drug effect than the spectral edge frequency (SE95
), we chose to use the latter for consistency with previous work investigating the influence of age on opioid‐induced EEG changes. Based on the spectral edge as a measure of opioid drug effect, we reported an age‐related increase in the half‐life ke0
for remifentanil. This slower equilibration between blood and brain concentrations in the elderly will offset their higher initial blood concentrations. Consequently, peak effect site concentrations after a bolus dose will be similar in the elderly and in the young, although the time of the peak will be relatively delayed. Thus, based on the EEG model, there are no pharmacokinetic grounds for recommending reduced bolus doses in the elderly, and it is for pharmacodynamic reasons (the 50% reduction in EC50
in the elderly) that we recommend halving bolus doses in the elderly.
The maintenance infusion rate is calculated as the product of the target concentration and the clearance of the drug. Again, based on the EEG findings of a 50% decrease in EC50 in elderly participants, we can anticipate that there will be at least a 50% decrease in maintenance infusion rates in elderly persons. However, the clearance in the elderly is about 66% of clearance in younger persons. This further decreases the infusion rate required to maintain any given drug effect. The net effect of the change in remifentanil potency and clearance in elderly patients is that the infusion rate to maintain any level of drug effect is about one third the rate in a younger person (66% x 50% = 33%).
We did not use the same weighted measures of performance to evaluate the pharmacokinetic and pharmacodynamic data. The pharmacodynamic data spanned a range from approximately 5 to 20 Hz. In contrast to the pharmacokinetic data (Figure 2
) this range is less than one order of magnitude. Within this range the variance of the residual intraindividual error was well described by an unscaled error model. Because we used an unscaled error model for the objective function, we believe it more appropriate to use an unscaled measure for the residual error. However, if we had used the same weighted measures of performance, we would have found that the MDAWR and the MAWR of the pharmacodynamic data were only slightly greater than for the pharmacokinetic data. Specifically, the MDAWR and MAWR were 15.3% and 22.8%, respectively, for the complex pharmacokinetic model (Table 3
) and were 17.2% and 24.7%, respectively, for the complex pharmacodynamic model.
Our finding of an increased sensitivity in the elderly using the EEG model is consistent with the earlier work of Scott and Stanski 
for fentanyl and alfentanil. Although the EEG is not, in and of itself, a clinically useful end point, it offers several experimental advantages over clinical measures of drug effect. The EEG is a nearly instantaneous measure of drug effect, a continuous measure of drug effect, and an objective and reproducible measure of drug effect. These benefits of using the EEG would be irrelevant if the EEG told us nothing about the time course of clinically important opioid drug effects or about the relative potency of the opioids. The EEG predicts a more rapid and evanescent effect from a bolus of remifentanil and alfentanil compared with fentanyl or sufentanil. This correlates with the observed time course of effect of these drugs in the operating room for such effects as ventilatory depression, sedation, and analgesia. Shafer and Varvel 
reviewed the pharmacodynamics of fentanyl, alfentanil, and sufentanil and observed that the EEG predicted relative potencies among these drugs that were entirely consistent with the predicted EC50
or therapeutic ranges for clinical measures of drug effect. Hill and associates 
examined analgesic response, side effects, and evoked potentials during steady‐state opioid infusions and found a robust relation between plasma drug concentration and analgesic, ventilatory, and subjective‐effect magnitudes for each opioid in the study. We are aware of no data that suggest that the relative potency for opioids as measured by the EEG differs from the relative potency based on clinical measures of drug effect. The opioid studies that did not observe changes in brain sensitivity with age, those of Lemmens and coworkers 
and Ausems and colleagues, 
used low‐resolution pharmacodynamic end points and examined fewer participants. Consequently, their studies may have lacked the resolution and the statistical power to detect the changes we report. However, the EEG is not clinical anesthesia, and it is possible that changes in potency reflected in the EEG may not reflect changes in potency for some clinical end points.
Although statistically significant, the addition of age and LBM as covariates resulted in a relatively small improvement in the prediction of EEG effect when applied prospectively. This may have been due in part to the study design. Relatively high infusion rates were chosen to achieve maximum EEG effect in all participants for the purpose of developing the full pharmacokinetic and pharmacodynamic model. Thus the EEG rapidly changed from baseline to high amplitude, low‐frequency delta waves in all participants. As neither E0
were related to subject covariates, it is not surprising that the simple and complex pharmacodynamic models performed similarly during the infusion. However, even with this study design, covariate‐related differences in pharmacokinetics and pharmacodynamics should allow a more accurate prediction of the EEG effect during the postinfusion period. In the group used to derive the models (Figure 7
) and the group used for the prospective test, the pattern of the residuals for the pharmacodynamic data improved during the postinfusion period, although this was relatively small. Thus, even when statistically significant relations between model parameters and participant covariates are detected, further analysis is required to evaluate the importance of these relations in the presence of known interindividual variability in model parameters.
We used NONMEM and GAM to develop simple and complex population pharmacokinetic and pharmacodynamic models, and we assessed these models prospectively. We found (1) an effect of age on the pharmacokinetic and pharmacodynamic parameters for remifentanil; (2) an effect of LBM on the pharmacokinetic parameters; and (3) no effect of gender on any pharmacokinetic or pharmacodynamic parameter. In our companion article, we will use the estimates of interindividual variability in pharmacokinetic and pharmacodynamic parameters provided by NONMEM to perform computer simulations that examine the influence of age, LBM, and intersubject variability on remifentanil dosing guidelines.
*Beal SL, Sheiner LB: NONMEM User's Guide. San Francisco, University of California San Francisco, 1979. The pharmacokinetic and pharmacodynamic models for the NONMEM and GAM analyses are available via the WWW at URL: http://pkpd.icon.palo‐alto.med.va.‐gov.
**S‐PLUS for Windows, Version 3.2, copyright 1988, 1994 Mathsoft, Inc.
***James WPT: Research on obesity. London, Her Majesty's Stationary Office, 1976.
1. Lemmens HJ, Burm AG, Hennis PJ, Gladines MP, Bovill JG: Influence of age on the pharmacokinetics of alfentanil. Gender dependence. Clin Pharmacokinet 1990; 19:416-22.
2. Scott JC, Stanski DR: Decreased fentanyl and alfentanil dose requirements with age. A simultaneous pharmacokinetic and pharmacodynamic evaluation. J Pharmacol Exp Ther 1987; 240:159-66.
3. Lemmens HJ, Bovill JG, Hennis PJ, Burm AG: Age has no effect on the pharmacodynamics of alfentanil. Anesth Analg 1988; 67:956-60.
4. Lemmens HJ, Burm AG, Bovill JG, Hennis PJ: Pharmacodynamics of alfentanil as a supplement to nitrous oxide anaesthesia in the elderly patient. Br J Anaesth 1988; 61:173-9.
5. Egan TD, Lemmens HJ, Fiset P, Hermann DJ, Muir KT, Stanski DR, Shafer SL: The pharmacokinetics of the new short-acting opioid remifentanil (GI87084B) in healthy adult male volunteers. Anesthesiology 1993; 79:881-92.
6. Egan TD, Minto CF, Hermann DJ, Barr J, Muir KT, Shafer SL: Remifentanil vs alfentanil: Comparative pharmacokinetics and pharmacodynamics in healthy adult male volunteers. Anesthesiology 1996; 84:821-33.
7. Grosse CM, Davis IM, Arrendale RF, Jersey J, Amin J: Determination of remifentanil in Human blood by liquid-liquid extraction and capillary GC-HRMS-SIM using a deuterated internal standard. J Pharm Biomed Anal 1994; 12:195-203.
8. Jersey JA, Guyan S, Abbey L, Grosse CM, Davis IM: Determination of Alfentanil in whole blood by GC/MS using liquid-liquid and solid phase extractions. Pharm Res 1993; 10(S):1180.
9. Mandema JW, Verotta D, Sheiner LB: Building population pharmacokinetic-pharmacodynamic models. I. Models for covariate effects. J Pharmacokinet Biopharm 1992; 20:511-28.
10. Sheiner LB, Beal S, Rosenberg B, Marathe VV: Forecasting individual pharmacokinetics. Clin Pharmacol Ther 1979; 26:294-305.
11. Sheiner LB, Beal SL: Bayesian individualization of pharmacokinetics: Simple implementation and comparison with non-Bayesian methods. J Pharm Sci 1982; 71:1344-8.
12. DuBois D, DuBois EF: Clinical calorimetry: X. A formula to estimate the approximate surface area if the height and weight be known. Arch Intern Med 1916; 17:863-71.
13. Hastie TJ: Generalized additive models, Statistical Models. Edited by Chambers JM, Hastie TJ. London, Chapman & Hall, 1993, pp 249-307.
14. Kataria BK, Ved SA, Nicodemus HF, Hoy GR, Lea D, Dubois MY, Mandema JW, Shafer SL: The pharmacokinetics of propofol in children using three different data analysis approaches. Anesthesiology 1994; 80:104-22.
15. Billard V, Gambus PL, Barr J, Minto CF, Corash L, Tessman JT, Stickney JL, Shafer SL: The pharmacokinetics of 8-methoxypsoralen following intravenous administration in humans. Br J Clin Pharm 1995; 40:347-60.
16. Scott JC, Cooke JE, Stanski DR: Electroencephalographic quantitation of opioid effect: Comparative pharmacodynamics of fentanyl and sufentanil. Anesthesiology 1991; 74:34-42.
17. Varvel JR, Donoho DL, Shafer SL: Measuring the predictive performance of computer-controlled infusion pumps. J Pharmacokinet Biopharm 1992; 20:63-94.
18. Minto CF, Schnider TW, Shafer SL: Pharmacokinetics and pharmacodynamics of remifentanil. II. Model application. Anesthesiology 1997; 86:24-33.
19. Helmers H, van Peer A, Woestenborghs R, Noorduin H, Heykants J: Alfentanil kinetics in the elderly. Clin Pharmacol Ther 1984; 36:239-43.
20. Hudson RJ, Thomson IR, Burgess PM, Rosenbloom M: Alfentanil pharmacokinetics in patients undergoing abdominal aortic surgery. Can J Anaesth 1991; 38:61-7.
21. Matteo RS, Schwartz AE, Ornstein E, Young WL, Chang WJ: Pharmacokinetics of sufentanil in the elderly surgical patient. Can J Anaesth 1990; 37:852-6.
22. Lehmann KA, Sipakis K, Gasparini R, van Peer A: Pharmacokinetics of sufentanil in general surgical patients under different conditions of anaesthesia. Acta Anaesthesiol Scand 1993; 37:176-80.
23. Helmers JH, van Leeuween L, Zuurmond WW: Sufentanil pharmacokinetics in young adult and elderly surgical patients. Eur J Anaesthesiol 1994; 11:181-5.
24. Ausems ME, Hug CC, Stanski DR, Burm AGL: Plasma concentrations of alfentanil required to supplement nitrous oxide anesthesia for general surgery. Anesthesiology 1986; 65:362-73.
25. Maitre PO, Buhrer M, Thomson D, Stanski DR: A three-step approach combining Bayesian regression and NONMEM population analysis: Application to midazolam. J Pharmacokinet Biopharm 1991; 19:377-84.
26. Westmoreland CL, Hoke JF, Sebel PS, Hug CC, Muir KT: Pharmacokinetics of remifentanil (GI87084B) and its major metabolite (GI90291) in patients undergoing elective inpatient surgery. Anesthesiology 1993; 79:893-903.
27. Gregg KM, Varvel JR, Shafer SL: Application of semilinear canonical correlation to the measurement of opioid drug effect. J Pharmacokinet Biopharm 1992; 20:611-35.
28. Gambus PL, Gregg KM, Shafer SL: Validation of the alfentanil canonical univariate parameter as a measure of opioid effect on the electroencephalogram. Anesthesiology 1995; 83:747-56.
29. Shafer SL, Varvel JR: Pharmacokinetics, pharmacodynamics, and rational opioid selection. Anesthesiology 1991; 74:53-63.
30. Hill HF, Chapman CR, Saeger LS, Bjurstrom R, Walter MH, Schoene RB, Kippes M: Steady-state infusions of opioids in human. II. Concentration-effect relationships and therapeutic margins. Pain 1990; 43:69-79.
© 1997 American Society of Anesthesiologists, Inc.