Free
Perioperative Medicine  |   June 2017
An Allometric Model of Remifentanil Pharmacokinetics and Pharmacodynamics
Author Notes
  • Department of Anesthesiology, University Medical Center Groningen, University of Groningen, The Netherlands (D.J.E., J.H.P., H.V., A.R.A., M.M.R.F.S.); Department of Anesthesiology, Leiden University Medical Center, Leiden, The Netherlands (E.O., J.V.); and Department of Anesthesia, Ghent University, Gent, Belgium (M.M.R.F.S.).
  • Submitted for publication October 29, 2015. Accepted for publication August 9, 2016.
    Submitted for publication October 29, 2015. Accepted for publication August 9, 2016.×
  • Corresponding articles on pages 993 and 1019.
    Corresponding articles on pages 993 and 1019.×
  • Supplemental Digital Content is available for this article. Direct URL citations appear in the printed text and are available in both the HTML and PDF versions of this article. Links to the digital files are provided in the HTML text of this article on the Journal’s Web site (www.anesthesiology.org).
    Supplemental Digital Content is available for this article. Direct URL citations appear in the printed text and are available in both the HTML and PDF versions of this article. Links to the digital files are provided in the HTML text of this article on the Journal’s Web site (www.anesthesiology.org).×
  • Address correspondence to Dr. Eleveld: Department of Anesthesiology, University Medical Center Groningen, Hanzeplein 1, 9713 GZ, Groningen, The Netherlands. d.j.eleveld@umcg.nl. Information on purchasing reprints may be found at www.anesthesiology.org or on the masthead page at the beginning of this issue. Anesthesiology’s articles are made freely accessible to all readers, for personal use only, 6 months from the cover date of the issue.
Article Information
Perioperative Medicine / Clinical Science / Pain Medicine / Pharmacology
Perioperative Medicine   |   June 2017
An Allometric Model of Remifentanil Pharmacokinetics and Pharmacodynamics
Anesthesiology 6 2017, Vol.126, 1005-1018. doi:10.1097/ALN.0000000000001634
Anesthesiology 6 2017, Vol.126, 1005-1018. doi:10.1097/ALN.0000000000001634
Abstract

Background: Pharmacokinetic and pharmacodynamic models are used to predict and explore drug infusion schemes and their resulting concentration profiles for clinical application. Our aim was to develop a pharmacokinetic-pharmacodynamic model for remifentanil that is accurate in patients with a wide range of age and weight.

Methods: Remifentanil pharmacokinetic data were obtained from three previously published studies of adults and children, one of which also contained pharmacodynamic data from adults. NONMEM was used to estimate allometrically scaled compartmental pharmacokinetic and pharmacodynamic models. Weight, age, height, sex, and body mass index were explored as covariates. Predictive performance was measured across young children, children, young adults, middle-aged, and elderly.

Results: Overall, 2,634 remifentanil arterial concentration and 3,989 spectral-edge frequency observations from 131 individuals (55 male, 76 female) were analyzed. Age range was 5 days to 85 yr, weight range was 2.5 to 106 kg, and height range was 49 to 193 cm. The final pharmacokinetic model uses age, weight, and sex as covariates. Parameter estimates for a 35-yr-old, 70-kg male (reference individual) are: V1, 5.81 l; V2, 8.82 l; V3, 5.03 l; CL, 2.58 l/min; Q2, 1.72 l/min; and Q3, 0.124 l/min. Parameters mostly increased with fat-free mass and decreased with age. The pharmacodynamic model effect compartment rate constant (ke0) was 1.09 per minute (reference individual), which decreased with age.

Conclusions: We developed a pharmacokinetic-pharmacodynamic model to predict remifentanil concentration and effect for a wide range of patient ages and weights. Performance exceeded the Minto model over a wide age and weight range.

What We Already Know about This Topic
  • A widely used remifentanil pharmacokinetic-pharmacodynamic model was developed from a study of healthy adult volunteers

  • Allometry, which addresses the relationship of body size to shape, anatomy, physiology, and behavior, might facilitate development of a model that could be useful in both children and adults

What This Article Tells Us That Is New
  • A general-purpose remifentanil pharmacokinetic-pharmacodynamic model was developed using pharmacokinetic data from studies of adults and children and pharmacodynamic data from an adult study

  • Model parameters were influenced by the patient covariates fat-free mass, weight, age, and sex

  • The predictive performance of the model was in a clinically acceptable range for all subgroups considered and was better than that of a widely used model, particularly in young children and children

PHARMACOKINETIC (PK) and pharmacodynamic (PD) models are used to predict drug concentration and effect profiles from a given drug administration scheme and to predict the drug administration scheme necessary to achieve a desired drug effect. They are useful for intraoperative advisory displays1  to provide clinicians with feedback on expected drug concentrations and effects. Target-controlled infusion (TCI) systems use PK-PD models to calculate the infusion rates necessary to achieve desired drug concentrations or effects.
For remifentanil, the PK-PD model developed by Minto et al.2  is widely used. It was derived from a study of 65 healthy adult volunteers aged from 20 to 85 yr with body weights from 45 to 106 kg. Others have studied remifentanil PK in children3  and adults.4  Combining these data might lead to a more general applicable PK-PD model for remifentanil. A similar approach was taken by Eleveld et al.,5  who developed a general-purpose PK model for propofol by combining data from multiple studies from young children, children, adults, the obese, and the elderly. The resulting model performed better than specialized models and has been prospectively validated.6–8 
Because of the quite different sizes of children and adults, concepts from allometry are likely to be beneficial to model development. Allometry concerns the relationship of body size to shape, anatomy, physiology, and behavior. While allometry has been applied in the biologic sciences dating back to the 1930s, with Kleiber’s law,9  more recently it has been systematically applied to PK models,10  and it has been used for remifentanil.11  Typically, it involves the initial hypothesis that volumes scale linearly with body size (usually but not always weight) and clearances to body size to the ¾ power.
Our aim was to combine several remifentanil PK and PD data sets to develop a combined PK-PD model. We hypothesized that this may lead to improved predictive performance compared to previously published models.
Materials and Methods
We obtained remifentanil PK and PD data from three studies: the aforementioned study in adults by Minto et al.,2  a study in children by Ross et al.,3  and a study in adults with TCI by Mertens et al.4  For the studies by Minto et al. and Ross et al., the original data were available via the Open TCI Initiative web site (http://www.opentci.org). For the study by Mertens et al., the remifentanil TCI infusion targets, PK model, patient characteristics, and target concentration data were provided to us by the authors. These data allow the infusion profile in each individual to be back-calculated and validated with the predicted concentrations that were recorded during performance of the study. Table 1 provides some summarized details of the data sets considered. All of the data sets have been published in scientific journals, and the necessary ethical committee approval was obtained before performance of the studies.
Table 1.
Details of the Component Data Sets
Details of the Component Data Sets×
Details of the Component Data Sets
Table 1.
Details of the Component Data Sets
Details of the Component Data Sets×
×
In the data set, height information is missing for some individuals, and for others the recorded height seems unreasonable. For these individuals, their height was assumed to be the average height of other individuals whose weight was within 10% of their weight. Model estimation and evaluation were performed using NONMEM version 7.3 (Icon Development Solutions, USA) using the first order conditional estimation method with interaction. Calculations of predictive performance were performed using the R language,12  version 2.14.1. Simulations of TCI using the models were performed using PK-PD Tools for Excel© (www.pkpdtools.com).
PK and PD Analysis
We modeled the time course of remifentanil plasma concentration using a three-compartment PK model with volumes V1, V2, V3, elimination clearance CL, and intercompartmental clearances Q2 and Q3. Our initial model scaled all parameters linearly with total body weight (TBW). We also tested allometric scaling where volumes are scaled linearly with body size and clearances to a power exponent of 0.75. These are sometimes referred to as theoretical values for scaling exponents.10  We also considered the fat-free mass (FFM) predictor described by Al-Sallami et al.13  as a body-size descriptor that uses a maturation model to extend the Janmahasatian et al.14  adult FFM predictor to include children. The Al-Sallami FFM predictor uses TBW (kg), age (yr), body mass index (BMI), and sex:
Residual error in PK observations was assumed to be log-normally distributed using the “transform-both-sides” approach (with the method NONMEM uses to model intrasubject variability, true lognormal distributions of residual errors can only be modeled by comparing the log of the predicted concentration with the log of the measured concentration using a simple additive error model), and a separate residual error was estimated for each data set.
We modeled PD measures using a sigmoidal Emax model driven by an effect compartment concentration (Ce) connected to the plasma compartment by a first-order rate constant (ke0). Model parameters were assumed to be log-normally distributed across the population. Residual error in PD observations was assumed additive and normally distributed. The equation of the PD model was as follows,
where C and Ce are the concentrations in the central (V1) and effect compartments, E0 is the baseline PD measure when no drug is present, Emax is the maximum possible drug effect, Ce50 is the Ce associated with 50% of the maximum effect, γ is the steepness of the concentration versus response relation, and ε represents additive residual error to the PD observations. For PD model estimation, the individual predicted plasma concentrations from the final PK model were used as the driving force for the effect compartment; this is known as the sequential method,15,16  also known as the individual post hoc prediction method.
Model parameters were calculated relative to a reference individual,17  a 70-kg, 35-yr, 170-cm male. Age, weight, height, sex, and BMI were considered as potential covariates for model parameters. Potential covariate relationships were identified by examination of individual variability (η) values obtained from NONMEM estimation and tested for inclusion in the model. Parameters and covariate relationships were added and removed from the model during hierarchical model building to obtain a good model fit using the corrected Akaike information criterion (AIC). We required a decrease in AIC of at least 10 for adding parameters to the model and allowed an increase in AIC of up to 4 when removing parameters. We calculated simple measures of the ability of each pharmacokinetic and pharmacodynamic model to predict the observations. We did not consider models that show clear degradation of the ability of models to predict the observations, even if such models had improved NONMEM objective functions. Details of the measures are described in the following section. Overall, our approach corresponds to requiring strong evidence for adding parameters to the model and allowing moderate degradation in model fit for simplifications of the model. These criteria would generally be considered rather conservative with a tendency to result in a (comparatively) simple final model.
Uncertainty in estimated model parameters was evaluated by estimating the upper and lower 95% confidence limits by spline interpolation of the likelihood profiles. We determined what increase/decrease in each parameter is required to increase NONMEM objective function by 3.84. For estimates of logarithmic interindividual variability, we report the estimated variance and also the coefficient of variation using the method described by Elassaiss-Schaap and Heisterkamp,*
where ω2 is the estimated parameter population variance.
Quantifying PK and PD Predictive Performance
To quantify the PK predictive performance for an observation, we calculated the performance error18  (PEPK) and absolute performance error (APEPK) as follows.
For PD measures, only the Minto study data provided PD observations that were the spectral-edge frequency (SEF) of the electroencephalograph.19  For these, we used prediction error calculations appropriate for additive error models.
For these PK and PD performance error measures, the median values are reported: the median PEPK (MdPEPK), median APEPK (MdAPEPK), median PEPD (MdPEPD), and median APEPD (MdAPEPD). The MdPE measures indicate bias, and the MdAPE measures indicate precision.
During model development, the model predictive performance was evaluated using predictions from repeated twofold cross-validation. With this technique, before analysis, individuals with observation data are randomly partitioned into equally sized groups D1 and D2. First, the model parameters are estimated from the data of D1. The parameters are held constant and then used to predict the observations of group D2 using population predictions based on the measured covariates (η fixed to 0). The process is then repeated after exchanging D1 and D2. The predictions for D1 and D2 are combined to obtain a complete set of predictions. Thus predictive performance is only evaluated for observations that were not used for model estimation, i.e., out-of-sample predictions. In an effort to reduce Monte-Carlo variability in performance estimates caused by random partitioning, twofold cross-validation was repeated 10 times with different random partitions of individuals, and all predictions were combined and used to evaluate model predictive performance.
If observations in a data set are not balanced with respect to relevant subgroups, then the model that strictly minimizes AIC may not result in balanced performance for all subgroups of the data. This can occur when model modifications improve model fit for well-represented subgroups to the disproportionate detriment of less well-represented but clinically important subgroups. To address this, we split the individuals studied into five subgroups: young children (<3 yr), children (3 to <18 yr), young adults (18 to <40 yr), middle-aged (40 to <65 yr), and elderly (≥65 yr). The overall measure of predictive performance was the average MdAPEPK and MdAPEPD across these subgroups. The subgroups were chosen to approximately parallel PK models in the literature.
Comparison to Existing Models
The performance of the final model was compared to models available in the literature. For the models from the literature, no cross-validation was performed, i.e., in-sample predictions. The models considered were the Minto model,2  the Egan model,20  and the La Colla model.21  We also considered the Rigby-Jones model,11  which is an allometric PK model developed from data from children.
Results
The analyzed data set contains 2,634 remifentanil concentration observations from three studies and 3,989 SEF observations from one study. Overall, data from 131 individuals (55 male, 76 female) were studied. The age range was from 5 days to 85 yr, weight range was from 2.5 to 106 kg, and height range was from 49 to 193 cm. Figure 1 shows the distribution of age, weight, height, and BMI of the individuals. The distributions of individuals and observations for each of the subgroups considered are shown in table 2.
Table 2.
Composition of Subgroups Used for Intraoperative Predictive Performance Evaluation
Composition of Subgroups Used for Intraoperative Predictive Performance Evaluation×
Composition of Subgroups Used for Intraoperative Predictive Performance Evaluation
Table 2.
Composition of Subgroups Used for Intraoperative Predictive Performance Evaluation
Composition of Subgroups Used for Intraoperative Predictive Performance Evaluation×
×
Fig. 1.
Histograms for age, weight, height, body mass index (BMI), and weight versus age and height versus weight for the individuals studied.
Histograms for age, weight, height, body mass index (BMI), and weight versus age and height versus weight for the individuals studied.
Fig. 1.
Histograms for age, weight, height, body mass index (BMI), and weight versus age and height versus weight for the individuals studied.
×
For four individuals height information was missing. Two 12-yr-old children had recorded heights that seem unlikely given their weight and age, leading to BMI values of 38.8 and 106. Height values for these individuals were treated as missing, and the estimated height was used for FFM and BMI calculations.
PK Model Development
A graphical description of the PK model development process is shown in figure 2. The initial linear model (Model 1) has six basic PK parameters (V1, V2, V3, CL, Q2, and Q3), each scaled linearly with TBW and each with interindividual variance. There are also three residual error estimates, one for each component data set. Scaling all parameters by FFM predicted by the Al-Sallami model (Model 2) and the addition of allometric scaling (Model 3) both led to improved models, evidenced by lower AIC and improved predictive performance. The application of compartmental allometry5  to Q2 and Q3 (Model 4) further improved the model. In this approach, Q2 and Q3 scale to the power exponent 0.75 to the relative estimated size of V2 and V3. Here the size scaling factors of Q2 and Q3 were (V2/V2ref)0.75 and (V3/V3ref)0.75, respectively instead of the usual allometric approach of scaling with TBW, i.e., with (TBW/TBWref)0.75.
Fig. 2.
Hierarchical pharmacokinetic model building. AIC = Akaike information criteria; CL = remifentanil clearance; k = number of model parameters; MdAPE = median absolute prediction error for all observations; Avg MdAPE = median (cross-validation) absolute prediction error averaged over the subgroups young children, children, young adults, middle-aged, and elderly.
Hierarchical pharmacokinetic model building. AIC = Akaike information criteria; CL = remifentanil clearance; k = number of model parameters; MdAPE = median absolute prediction error for all observations; Avg MdAPE = median (cross-validation) absolute prediction error averaged over the subgroups young children, children, young adults, middle-aged, and elderly.
Fig. 2.
Hierarchical pharmacokinetic model building. AIC = Akaike information criteria; CL = remifentanil clearance; k = number of model parameters; MdAPE = median absolute prediction error for all observations; Avg MdAPE = median (cross-validation) absolute prediction error averaged over the subgroups young children, children, young adults, middle-aged, and elderly.
×
We tested an exponential age covariate relationship centered at 35 yr17  for all parameters (Model 5). This involves multiplying each parameter by a function of the form faging(x)=exp(x · (AGE - 35)) and estimating different values for x for each parameter. This resulted in a large decrease in AIC, but with a small (0.1%) decrease in predictive performance. We provisionally accepted this into the model and planned to consider simplification of the aging model later in development. For Model 5, there was systematic deviation in CL for small-sized individuals. The addition of an sigmoidal maturation function of the form fsigmoid(x,E50,λ = xλ/(xλ + E50λ), where x is post- menstrual-age (PMA), an approach advocated by Anderson and Holford,10,22  led to an improvement in AIC (Model 6) and predictive performance as measured by the MdAPEPK averaged across the patient subgroups also improved. Using TBW as x (Model 7) instead of PMA led to a greater improvement in AIC and predictive performance, and so we accepted this feature into the model. The slope parameter λ could not be reliably estimated and was empirically fixed to 2. This corresponds to a smooth gradual increase of CL in early life up to adult values achieved as maturation completes.
Clearance was found to be increased in females (Model 8); however, increasing V2, CL, and Q2 for females aged approximately 12 to 45 yr (Model 9) showed greater improvement in AIC and predictive performance. Estimating these 12- and 45-yr boundaries from the data (Model 10) led to a small improvement in objective function but lowered predictive performance, suggesting overfitting. The slope parameter for the boundaries could also not be reliably estimated and was empirically fixed to 6, which provides a smooth transition.
Model 9 showed significant deviation in V3 from the typical value, suggesting that the assumed allometric relationship with body size is inadequate. Adding an exponential weight correction to V3 (Model 11) led to an improved model. Further, the age relationships for V1, Q2, and Q3 were similar, as were those for CL and V2. Estimating shared age covariates for these parameters (Model 12) reduced the number of parameters in the model and improved AIC, and predictive performance was unchanged. Other model modifications were tried, but none improved both AIC and predictive performance. We did evaluate other size-scaling methods, but none resulted in a better final model. An exploration of these results is provided in the supplementary data (Supplemental Digital Content 1, http://links.lww.com/ALN/B421).
The complete PK data set and the full NONMEM model code of the final PK model can be found in Supplemental Digital Content 2 (http://links.lww.com/ALN/B422). The likelihood profiles (Supplemental Digital Content 3, http://links.lww.com/ALN/B423) and the relationships between η and the covariates (Supplemental Digital Content 4, http://links.lww.com/ALN/B424) and between η1 to η6 (Supplemental Digital Content 5, http://links.lww.com/ALN/B425) can be found in the supplementary data, as well as a Visual Basic Macro suitable for with PK-PD Tools for Excel© (www.pkpdtools.com). The summarized equations of the final model are as follows.
Constants with the subscript ref are calculated for the reference individual. Symbols V1ref, V2ref, V3ref, CLref, Q2ref, and Q3ref are the estimated compartmental volumes and clearances for the reference individual, Θ1 defines the weight at which 50% of maturation of CL is complete, Θ2 defines the change in V1, Q2, and Q3 with age, Θ3 defines the change in V2 and CL with age, Θ4 defines the change in V3 with age, Θ5 defines the change in V2, CL, and Q2 for females aged 12 to 45 yr, and Θ6 defines the deviation of V3 from theoretical allometric scaling. Symbols η1 to η6 represent random variables with variances denoted in table 3, where η5 and η6 represent the intersubject variability of Q2 and Q3 not captured by that of V2 and V3. The symbols AGE and TBW represent an individual’s age in yr and weight in kg, respectively. Symbol ε represents residual observation error in the log domain with a variance fixed to 1. RES is residual standard deviation estimated separately for each data set in table 1.
Table 3.
Estimated Model Parameters and Population Variances in the Final Pharmacokinetic Model
Estimated Model Parameters and Population Variances in the Final Pharmacokinetic Model×
Estimated Model Parameters and Population Variances in the Final Pharmacokinetic Model
Table 3.
Estimated Model Parameters and Population Variances in the Final Pharmacokinetic Model
Estimated Model Parameters and Population Variances in the Final Pharmacokinetic Model×
×
Figure 3 shows population and individual predictions versus time and observed remifentanil plasma concentrations for the entire data set. Figure 4 shows the individual estimated PK model parameters plotted against weight. Although V1, V2, CL, and Q2 show what may be described as classical behavior, increasing with weight and decreasing with age for adults, V3 and Q3 show considerably more variability and appear to decrease with increasing weight above about 20 kg. Figure 5 shows the individual estimated PK model parameters plotted against age. The decline in parameter values with advancing age is visible for most parameters.
Fig. 3.
Population and post hoc pharmacokinetic predictions for the current study versus time and observed remifentanil plasma concentration. The black line is a Loess smoother.
Population and post hoc pharmacokinetic predictions for the current study versus time and observed remifentanil plasma concentration. The black line is a Loess smoother.
Fig. 3.
Population and post hoc pharmacokinetic predictions for the current study versus time and observed remifentanil plasma concentration. The black line is a Loess smoother.
×
Fig. 4.
Post hoc estimated pharmacokinetic model parameters plotted versus weight. Closed triangles = males; open circles = females. CL = remifentanil clearance.
Post hoc estimated pharmacokinetic model parameters plotted versus weight. Closed triangles = males; open circles = females. CL = remifentanil clearance.
Fig. 4.
Post hoc estimated pharmacokinetic model parameters plotted versus weight. Closed triangles = males; open circles = females. CL = remifentanil clearance.
×
Fig. 5.
Post hoc estimated pharmacokinetic model parameters plotted versus age. Closed triangles = males; open circles = females.
Post hoc estimated pharmacokinetic model parameters plotted versus age. Closed triangles = males; open circles = females.
Fig. 5.
Post hoc estimated pharmacokinetic model parameters plotted versus age. Closed triangles = males; open circles = females.
×
PK Performance Evaluation
Table 4 shows the estimation of predictive performance of the final PK model compared to other remifentanil PK models from the literature. The final PK model performed best for all groups except for young adults, for whom the in-sample performance was slightly lower than the Minto model (MdAPE 15.8 vs. 15.3%). Notably, the out-of-sample performance of the final model was equal to or better than the in-sample performance of the Minto model for the other age groups. The performance of the final PK model appears to exceed that of the Minto model over a wide age and weight range.
Table 4.
Estimates Predictive Performance for Pharmacokinetic (PK) Models Considered
Estimates Predictive Performance for Pharmacokinetic (PK) Models Considered×
Estimates Predictive Performance for Pharmacokinetic (PK) Models Considered
Table 4.
Estimates Predictive Performance for Pharmacokinetic (PK) Models Considered
Estimates Predictive Performance for Pharmacokinetic (PK) Models Considered×
×
Interestingly, the performance of the allometrically scaled Rigby-Jones model in adults is comparable to other models despite the fact that the model was developed from data from children. These results suggest that the primary difference between adults and children for remifentanil PK is simply size, for which the allometric model compensates with the scaling exponents. The models developed without allometric scaling on all parameters, the Minto, Egan, and La Colla models, performed poorly for children and very poorly for young children. These models do not adequately compensate for the size of the individual.
PD Model Development
For PD model development, Ce50 and ke0 were assumed independent of size and log-normally distributed across the population. The initial PD model showed correlation in η for ke0 with age, and adding an exponential age covariate led to improvement in AIC and predictive performance. Testing an exponential age covariate to Ce50, mirroring the PD model structure found by Minto, led to a small improvement in AIC (−6.92) but decreased predictive performance and was omitted from the final model. Restricting the age covariate only for elderly individuals also did not improve model performance. The complete PD data set and the full NONMEM model code are in Supplemental Digital Content 2 (http://links.lww.com/ALN/B422). The likelihood profiles can be found in Supplemental Digital Content 6 (http://links.lww.com/ALN/B426). The summarized equations of the final models are as follows.
Symbols E0ref, Emaxref, Ce50ref, γref, ke0ref, and Θ1 represent estimated model parameters, and η1 to η5 represent variances; the estimated values are shown in table 5. Figure 6 shows population and post hoc PD predictions for the current study versus time and observed SEF. Figures 7 and 8 show post hoc estimated PD model parameters ke0 and Ce50 plotted versus weight and age, respectively. A clear decrease with age is visible for ke0. For Ce50, age seems to most strongly influence elderly individuals, although that characteristic did not achieve statistical significance.
Table 5.
Estimated Model Parameters and Population Variances in the Final Pharmacodynamic (PD) Model
Estimated Model Parameters and Population Variances in the Final Pharmacodynamic (PD) Model×
Estimated Model Parameters and Population Variances in the Final Pharmacodynamic (PD) Model
Table 5.
Estimated Model Parameters and Population Variances in the Final Pharmacodynamic (PD) Model
Estimated Model Parameters and Population Variances in the Final Pharmacodynamic (PD) Model×
×
Fig. 6.
Population and post hoc pharmacodynamic predictions for the current study versus time and observed spectral-edge frequency (SEF). The black line is a Loess smoother.
Population and post hoc pharmacodynamic predictions for the current study versus time and observed spectral-edge frequency (SEF). The black line is a Loess smoother.
Fig. 6.
Population and post hoc pharmacodynamic predictions for the current study versus time and observed spectral-edge frequency (SEF). The black line is a Loess smoother.
×
Fig. 7.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus weight. Closed triangles = males; open circles = females.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus weight. Closed triangles = males; open circles = females.
Fig. 7.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus weight. Closed triangles = males; open circles = females.
×
Fig. 8.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus age. Closed triangles = males; open circles = females.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus age. Closed triangles = males; open circles = females.
Fig. 8.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus age. Closed triangles = males; open circles = females.
×
Discussion
We estimated a three-compartment allometric PK model with an effect-site compartment and sigmoidal Emax PD (spectral edge) model for remifentanil. The parameters of the model were influenced by patient covariates FFM, weight, age, and sex. For all subgroups, performance of the final model is in the range considered clinically acceptable (MdPE less than 10 to 20% and MdAPE between 20 and 40%).23,24 
Our PK model performed better than the Minto model, which is incorporated into commercially available TCI pumps and anesthesia displays, predicting a more diverse population to better accuracy while requiring fewer estimated parameters. The out-of-sample predictive performance of our model was equal to or better than the in-sample performance of the Minto model for all subgroups except young adults, for whom the performance was quite similar. We focused on balanced performance across subgroups by requiring model development to improve the average MdAPE over five subgroups. The cost of such a balanced approach can be slightly reduced performance in most well-represented subgroups, but the benefit is a more robust model.
Influence of Maturation
We found decreased CL for the smallest and youngest individuals in the data set. Using TBW as a predictor for maturation provided a slightly better model fit than PMA, the approach advocated by Anderson and Holford,10,22  although the differences were not very large. The final model estimated that 50% of maturation was complete at a bodyweight of approximately 3.6 kg, which is close to average birth weight, and maturation is 95% complete at approximately 14 kg. We cannot exclude the possibility that the smallest individuals had conditions associated with reduced remifentanil clearance.
Predictive Performance in the Obese
We did not analyze data from obese individuals, and our model must be extrapolated for this group. If the difference between the obese and nonobese is primarily size, then our model might be expected to perform well for the obese, because our model incorporates a size correction. However, if obesity entails physiologic changes overshadowing those attributable to size, then our model may extrapolate poorly. In contrast to the Minto model, we do not the use of the James equation for lean body mass; thus we avoid its paradoxical behavior for obese individuals.25  Additional data are needed to determine the performance of our model in the obese.
Influence of Sex
We found that sex indirectly influences the PK of remifentanil through its relationship with FFM. The influence of sex on V2, CL, and Q2 seems to be more complex, and we found that increasing these for females between 12 and 45 yr old resulted in an improved model. Broadly, this age range corresponds to the reproductive period in females, although the underlying physiologic mechanisms are not clear. The age range was obtained by visual inspection, and the data were not sufficiently informative to estimate the limits or the steepness at its boundaries. The model structure used assumes a smooth transition at the limits of the age range.
Allometric Scaling
We found that the theoretical scaling exponents of 1 for volumes and 0.75 for clearances were suitable for V1, V2, CL, Q2, and Q3. Only V3 deviated from our initial assumption of linear scaling with body size, showing an overall decrease in V3 for larger body sizes, and this is visible in figure 4. It is not clear what mechanisms play a role here. Although the non-specific esterases involved in remifentanil elimination26  can be found in diverse body tissues, there may be structural differences in large body sizes that are not accounted for in our model.
It is interesting that the Rigby-Jones model was developed in children but extrapolates reasonably well to adults. Conversely, all of the PK models tested that did not scale all parameters using allometry extrapolated poorly to children and very poorly on young children.
A limitation of the PD data is that it is not informative for potential allometric relations for ke0, due to the availability of adult data only. Allometric theory suggests physiologic rates scale to the −0.25 power exponent (or 0.25 if expressed as a half-life27 ). Heart rates and breathing rates are well known to scale allometrically for a wide range of body sizes. If the same applies to keo (it is also a physiologic rate), then one would expect it to scale to body size to the −0.25 power exponent. Testing this in the final PD model was inconclusive, with an AIC of 0.90 higher and predictive performance unchanged. Estimated parameters were only slightly changed or were identical to three significant digits. So allometric scaling for ke0 is not supported by the data, but at the same time, neither is the size independence of keo in the final model. Informative PD data are needed to resolve this issue. This likely requires PD data from the extremes of the size spectrum.
Compartmental Allometry
As reported for propofol,5  the application of compartmental allometry to the scaling of Q2 and Q3 led to improved model fit and predictive performance. Based on allometric theory, Anderson and Holford10  have proposed that compartmental volumes scale linearly with size. The reverse should also be true, that size scales linearly with volume. The insight of compartmental allometry is to use the estimated volumes as estimates of size (up to a constant factor) for the peripheral compartments. Thus using the individual estimates of V2 and V3 as size descriptors for Q2 and Q3 can be seen as an attempt to scale Q2 and Q3 by the estimated size of the compartment in the particular individual. The current study provides evidence that this approach performs better than using TBW as a size descriptor for peripheral compartments.
Effect-site Targeting
Figure 9 shows the remifentanil bolus doses and infusion rates needed to achieve 4 ng/ml in the effect compartment for the PK and PD models and for the Minto PK-PD model. Compared to the Minto model, our final model predicts a smaller initial bolus dose for 35-yr-old individuals and a similar dose in 70 yr-old male individuals. After the initial bolus, dose rates are quite similar. In the elderly, the time to the target effect is delayed by about 1 min compared to younger adults.
Fig. 9.
Remifentanil cumulative doses and infusion rates required to achieve 4 ng/ml in the effect compartment for the final model and the Minto model for 35- and 70-yr-old, 70-kg, 170-cm males.
Remifentanil cumulative doses and infusion rates required to achieve 4 ng/ml in the effect compartment for the final model and the Minto model for 35- and 70-yr-old, 70-kg, 170-cm males.
Fig. 9.
Remifentanil cumulative doses and infusion rates required to achieve 4 ng/ml in the effect compartment for the final model and the Minto model for 35- and 70-yr-old, 70-kg, 170-cm males.
×
Limitations of the Study
Prospective evaluation of the model was not performed. If different subgroups or measures of performance had been used, then a different model structure might be considered optimal. It could be argued that SEF has a tenuous relationship with other measures of remifentanil drug effect such as heart rate, arterial pressure, respiratory rate, tidal volume, muscle rigidity, and analgesia. One source of these doubts is that the Ce50 for SEF is 12.7 ng/ml, which is higher than TCI targets in routine clinical applications, which are typically in the range 3 to 8 ng/ml. Also, PD measurements were only available from adults from a single data set, so the PD parameters and the covariate relationships are likely to be less accurately identified than those from the PK. As such, our estimation of PD parameters such as ke0 and predictions of effect-site concentrations in children are extrapolations based on adult data.
The SEF data used to model the PD originate from the Minto model, which has been successfully applied in clinical practice for years, for effect-compartment-targeted remifentanil TCI in adults. As such, our model supports a wider covariate range for PK and a similar range for PD compared to the Minto model. We believe that clinicians using effect-compartmental-controlled TCI will benefit more from an integrated PK-PD model instead of adding a PD model to our PK model via indirect methods such as time to peak effect modeling.28 
Conclusion
We combined remifentanil PK data from adults and children and estimated a single PK-PD model, which shows good predictive performance for all subgroups considered. The safety and applicability of our final model needs to be evaluated.
Acknowledgments
This work builds on the efforts of many contributors. We used the PK and PD data obtained much from the Open TCI Initiative web site (http://www.opentci.org) collected through the efforts of Steven L. Shafer, M.D. (Department of Anesthesiology, Perioperative and Pain Medicine, Stanford, California), Charles F. Minto, M.B.Ch.B., F.A.N.Z.C.A., Ph.D. (Department of Anaesthesia, North Shore Private Hospital, Hunters Hill Private Hospital, Westmead Private Hospital, and Sydney Adventist Hospital, Sydney, Australia), Thomas W. Schnider, M.D. (Department for Anesthesiology, Intensive, Rescue and Pain Medicine, Kantonsspital St. Gallen, St. Gallen, Switzerland), and the numerous researchers who contributed data sets. The approach of open sharing of data has enabled us to focus on model building and leverage the immense work of those who planned the studies and collected the data sets. In turn, we share the complete data set analyzed along with NONMEM control files showing our final PK and PD models in the hope that other researchers may improve or extend on our results. We extend deep appreciation to all of the researchers, clinicians, patients, healthy volunteers, and others who directly and indirectly contributed.
Research Support
Support was provided solely from institutional and/or departmental sources.
Competing Interests
The authors declare no competing interests.
References
Struys, MM, Sahinovic, M, Lichtenbelt, BJ, Vereecke, HE, Absalom, AR . Optimizing intravenous drug administration by applying pharmacokinetic/pharmacodynamic concepts. Br J Anaesth 2011; 107:38–47 [Article] [PubMed]
Minto, CF, Schnider, TW, Egan, TD, Youngs, E, Lemmens, HJ, Gambus, PL, Billard, V, Hoke, JF, Moore, KH, Hermann, DJ, Muir, KT, Mandema, JW, Shafer, SL . Influence of age and gender on the pharmacokinetics and pharmacodynamics of remifentanil: I. Model development. Anesthesiology 1997; 86:10–23 [Article] [PubMed]
Ross, AK, Davis, PJ, Dear Gd, GL, Ginsberg, B, McGowan, FX, Stiller, RD, Henson, LG, Huffman, C, Muir, KT . Pharmacokinetics of remifentanil in anesthetized pediatric patients undergoing elective surgery or diagnostic procedures. Anesth Analg 2001; 93:1393–401, table of contents [Article] [PubMed]
Mertens, MJ, Olofsen, E, Engbers, FH, Burm, AG, Bovill, JG, Vuyk, J . Propofol reduces perioperative remifentanil requirements in a synergistic manner: Response surface modeling of perioperative remifentanil-propofol interactions. Anesthesiology 2003; 99:347–59 [Article] [PubMed]
Eleveld, DJ, Proost, JH, Cortínez, LI, Absalom, AR, Struys, MM . A general purpose pharmacokinetic model for propofol. Anesth Analg 2014; 118:1221–37 [Article] [PubMed]
Przybyłowski, K, Tyczka, J, Szczesny, D, Bienert, A, Wiczling, P, Kut, K, Plenzler, E, Kaliszan, R, Grześkowiak, E . Pharmacokinetics and pharmacodynamics of propofol in cancer patients undergoing major lung surgery. J Pharmacokinet Pharmacodyn 2015; 42:111–22 [Article] [PubMed]
Vereecke, HE, Eleveld, DJ, Colin, P, Struys, MM . Performance of the Eleveld pharmacokinetic model to titrate propofol in an obese Japanese patient population. Eur J Anaesthesiol 2016; 33:58–9 [Article] [PubMed]
Cortínez, LI, De la Fuente, N, Eleveld, DJ, Oliveros, A, Crovari, F, Sepulveda, P, Ibacache, M, Solari, S . Performance of propofol target-controlled infusion models in the obese: Pharmacokinetic and pharmacodynamic analysis. Anesth Analg 2014; 119:302–10 [Article] [PubMed]
Kleiber, M . Body size and metabolism. Hilgardia 1932; 6: 315–353. [Article]
Anderson, BJ, Holford, NH . Mechanism-based concepts of size and maturity in pharmacokinetics. Annu Rev Pharmacol Toxicol 2008; 48:303–32 [Article] [PubMed]
Rigby-Jones, AE, Priston, MJ, Sneyd, JR, McCabe, AP, Davis, GI, Tooley, MA, Thorne, GC, Wolf, AR . Remifentanil-midazolam sedation for paediatric patients receiving mechanical ventilation after cardiac surgery. Br J Anaesth 2007; 99:252–61 [Article] [PubMed]
R Development Core Team: R: A language and environment for statistical computing. 2008 Vienna, Austria, R Foundation for Statistical Computing, http://www.R-project.org.
Al-Sallami, HS, Goulding, A, Grant, A, Taylor, R, Holford, N, Duffull, SB . Prediction of fat-free mass in children. Clin Pharmacokinet 2015; 54:1169–78 [Article] [PubMed]
Janmahasatian, S, Duffull, SB, Ash, S, Ward, LC, Byrne, NM, Green, B . Quantification of lean bodyweight. Clin Pharmacokinet 2005; 44:1051–65 [Article] [PubMed]
Zhang, L, Beal, SL, Sheiner, LB . Simultaneous vs. sequential analysis for population PK/PD data I: Best-case performance. J Pharmacokinet Pharmacodyn 2003; 30:387–404 [Article] [PubMed]
Zhang, L, Beal, SL, Sheinerz, LB . Simultaneous vs. sequential analysis for population PK/PD data II: Robustness of methods. J Pharmacokinet Pharmacodyn 2003; 30:405–16 [Article] [PubMed]
Holford, NH . A size standard for pharmacokinetics. Clin Pharmacokinet 1996; 30:329–32 [Article] [PubMed]
Varvel, JR, Donoho, DL, Shafer, SL . Measuring the predictive performance of computer-controlled infusion pumps. J Pharmacokinet Biopharm 1992; 20:63–94 [Article] [PubMed]
Scott, JC, Cooke, JE, Stanski, DR . Electroencephalographic quantitation of opioid effect: Comparative pharmacodynamics of fentanyl and sufentanil. Anesthesiology 1991; 74:34–42 [Article] [PubMed]
Egan, TD, Huizinga, B, Gupta, SK, Jaarsma, RL, Sperry, RJ, Yee, JB, Muir, KT . Remifentanil pharmacokinetics in obese versus lean patients. Anesthesiology 1998; 89:562–73 [Article] [PubMed]
La Colla, L, Albertin, A, La Colla, G, Porta, A, Aldegheri, G, Di Candia, D, Gigli, F . Predictive performance of the ‘Minto’ remifentanil pharmacokinetic parameter set in morbidly obese patients ensuing from a new method for calculating lean body mass. Clin Pharmacokinet 2010; 49:131–9 [Article] [PubMed]
Anderson, BJ, Holford, NH . Mechanistic basis of using body size and maturation to predict clearance in humans. Drug Metab Pharmacokinet 2009; 24:25–36 [Article] [PubMed]
Schüttler, J, Kloos, S, Schwilden, H, Stoeckel, H . Total intravenous anaesthesia with propofol and alfentanil by computer-assisted infusion. Anaesthesia 1988; 43(suppl):2–7 [Article] [PubMed]
Glass, PS, Shafer, S, Reves, JG . Miller, RD . Intravenous drug delivery systems, Miller’s Anesthesia. Edited by 2005, pp Philadelphia, Pennsylvania, Elsevier ( Churchill Livingstone), 439–480.
Absalom, AR, Mani, V, De Smet, T, Struys, MM . Pharmacokinetic models for propofol: Defining and illuminating the devil in the detail. Br J Anaesth 2009; 103:26–37 [Article] [PubMed]
Davis, PJ, Stiller, RL, Wilson, AS, McGowan, FX, Egan, TD, Muir, KT . In vitro remifentanil metabolism: The effects of whole blood constituents and plasma butyrylcholinesterase. Anesth Analg 2002; 95:1305–7, table of contents [Article] [PubMed]
Anderson, BJ, Holford, NH, Woollard, GA, Chan, PL . Paracetamol plasma and cerebrospinal fluid pharmacokinetics in children. Br J Clin Pharmacol 1998; 46:237–43 [Article] [PubMed]
Minto, CF, Schnider, TW, Gregg, KM, Henthorn, TK, Shafer, SL . Using the time of maximum effect site concentration to combine pharmacokinetics and pharmacodynamics. Anesthesiology 2003; 99:324–33 [Article] [PubMed]
Fig. 1.
Histograms for age, weight, height, body mass index (BMI), and weight versus age and height versus weight for the individuals studied.
Histograms for age, weight, height, body mass index (BMI), and weight versus age and height versus weight for the individuals studied.
Fig. 1.
Histograms for age, weight, height, body mass index (BMI), and weight versus age and height versus weight for the individuals studied.
×
Fig. 2.
Hierarchical pharmacokinetic model building. AIC = Akaike information criteria; CL = remifentanil clearance; k = number of model parameters; MdAPE = median absolute prediction error for all observations; Avg MdAPE = median (cross-validation) absolute prediction error averaged over the subgroups young children, children, young adults, middle-aged, and elderly.
Hierarchical pharmacokinetic model building. AIC = Akaike information criteria; CL = remifentanil clearance; k = number of model parameters; MdAPE = median absolute prediction error for all observations; Avg MdAPE = median (cross-validation) absolute prediction error averaged over the subgroups young children, children, young adults, middle-aged, and elderly.
Fig. 2.
Hierarchical pharmacokinetic model building. AIC = Akaike information criteria; CL = remifentanil clearance; k = number of model parameters; MdAPE = median absolute prediction error for all observations; Avg MdAPE = median (cross-validation) absolute prediction error averaged over the subgroups young children, children, young adults, middle-aged, and elderly.
×
Fig. 3.
Population and post hoc pharmacokinetic predictions for the current study versus time and observed remifentanil plasma concentration. The black line is a Loess smoother.
Population and post hoc pharmacokinetic predictions for the current study versus time and observed remifentanil plasma concentration. The black line is a Loess smoother.
Fig. 3.
Population and post hoc pharmacokinetic predictions for the current study versus time and observed remifentanil plasma concentration. The black line is a Loess smoother.
×
Fig. 4.
Post hoc estimated pharmacokinetic model parameters plotted versus weight. Closed triangles = males; open circles = females. CL = remifentanil clearance.
Post hoc estimated pharmacokinetic model parameters plotted versus weight. Closed triangles = males; open circles = females. CL = remifentanil clearance.
Fig. 4.
Post hoc estimated pharmacokinetic model parameters plotted versus weight. Closed triangles = males; open circles = females. CL = remifentanil clearance.
×
Fig. 5.
Post hoc estimated pharmacokinetic model parameters plotted versus age. Closed triangles = males; open circles = females.
Post hoc estimated pharmacokinetic model parameters plotted versus age. Closed triangles = males; open circles = females.
Fig. 5.
Post hoc estimated pharmacokinetic model parameters plotted versus age. Closed triangles = males; open circles = females.
×
Fig. 6.
Population and post hoc pharmacodynamic predictions for the current study versus time and observed spectral-edge frequency (SEF). The black line is a Loess smoother.
Population and post hoc pharmacodynamic predictions for the current study versus time and observed spectral-edge frequency (SEF). The black line is a Loess smoother.
Fig. 6.
Population and post hoc pharmacodynamic predictions for the current study versus time and observed spectral-edge frequency (SEF). The black line is a Loess smoother.
×
Fig. 7.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus weight. Closed triangles = males; open circles = females.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus weight. Closed triangles = males; open circles = females.
Fig. 7.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus weight. Closed triangles = males; open circles = females.
×
Fig. 8.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus age. Closed triangles = males; open circles = females.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus age. Closed triangles = males; open circles = females.
Fig. 8.
Post hoc estimated pharmacodynamic model parameters ke0 and Ce50 plotted versus age. Closed triangles = males; open circles = females.
×
Fig. 9.
Remifentanil cumulative doses and infusion rates required to achieve 4 ng/ml in the effect compartment for the final model and the Minto model for 35- and 70-yr-old, 70-kg, 170-cm males.
Remifentanil cumulative doses and infusion rates required to achieve 4 ng/ml in the effect compartment for the final model and the Minto model for 35- and 70-yr-old, 70-kg, 170-cm males.
Fig. 9.
Remifentanil cumulative doses and infusion rates required to achieve 4 ng/ml in the effect compartment for the final model and the Minto model for 35- and 70-yr-old, 70-kg, 170-cm males.
×
Table 1.
Details of the Component Data Sets
Details of the Component Data Sets×
Details of the Component Data Sets
Table 1.
Details of the Component Data Sets
Details of the Component Data Sets×
×
Table 2.
Composition of Subgroups Used for Intraoperative Predictive Performance Evaluation
Composition of Subgroups Used for Intraoperative Predictive Performance Evaluation×
Composition of Subgroups Used for Intraoperative Predictive Performance Evaluation
Table 2.
Composition of Subgroups Used for Intraoperative Predictive Performance Evaluation
Composition of Subgroups Used for Intraoperative Predictive Performance Evaluation×
×
Table 3.
Estimated Model Parameters and Population Variances in the Final Pharmacokinetic Model
Estimated Model Parameters and Population Variances in the Final Pharmacokinetic Model×
Estimated Model Parameters and Population Variances in the Final Pharmacokinetic Model
Table 3.
Estimated Model Parameters and Population Variances in the Final Pharmacokinetic Model
Estimated Model Parameters and Population Variances in the Final Pharmacokinetic Model×
×
Table 4.
Estimates Predictive Performance for Pharmacokinetic (PK) Models Considered
Estimates Predictive Performance for Pharmacokinetic (PK) Models Considered×
Estimates Predictive Performance for Pharmacokinetic (PK) Models Considered
Table 4.
Estimates Predictive Performance for Pharmacokinetic (PK) Models Considered
Estimates Predictive Performance for Pharmacokinetic (PK) Models Considered×
×
Table 5.
Estimated Model Parameters and Population Variances in the Final Pharmacodynamic (PD) Model
Estimated Model Parameters and Population Variances in the Final Pharmacodynamic (PD) Model×
Estimated Model Parameters and Population Variances in the Final Pharmacodynamic (PD) Model
Table 5.
Estimated Model Parameters and Population Variances in the Final Pharmacodynamic (PD) Model
Estimated Model Parameters and Population Variances in the Final Pharmacodynamic (PD) Model×
×