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

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:

formula

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,

formula

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,*

formula

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.

formula

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.

formula

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
Fig. 1.

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.

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.

formula

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

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.

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.

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.

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

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.

formula

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
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. 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.

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.

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.

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

1.
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
2.
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
3.
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
4.
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
5.
Eleveld
DJ
,
Proost
JH
,
Cortínez
LI
,
Absalom
AR
,
Struys
MM
:
A general purpose pharmacokinetic model for propofol.
Anesth Analg
2014
;
118
:
1221
37
6.
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
7.
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
8.
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
9.
Kleiber
M
:
Body size and metabolism.
Hilgardia
1932
;
6
:
315
353
.
10.
Anderson
BJ
,
Holford
NH
:
Mechanism-based concepts of size and maturity in pharmacokinetics.
Annu Rev Pharmacol Toxicol
2008
;
48
:
303
32
11.
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
12.
R Development Core Team
:
R: A language and environment for statistical computing
.
Vienna, Austria
,
R Foundation for Statistical Computing
,
2008
13.
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
14.
Janmahasatian
S
,
Duffull
SB
,
Ash
S
,
Ward
LC
,
Byrne
NM
,
Green
B
:
Quantification of lean bodyweight.
Clin Pharmacokinet
2005
;
44
:
1051
65
15.
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
16.
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
17.
Holford
NH
:
A size standard for pharmacokinetics.
Clin Pharmacokinet
1996
;
30
:
329
32
18.
Varvel
JR
,
Donoho
DL
,
Shafer
SL
:
Measuring the predictive performance of computer-controlled infusion pumps.
J Pharmacokinet Biopharm
1992
;
20
:
63
94
19.
Scott
JC
,
Cooke
JE
,
Stanski
DR
:
Electroencephalographic quantitation of opioid effect: Comparative pharmacodynamics of fentanyl and sufentanil.
Anesthesiology
1991
;
74
:
34
42
20.
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
21.
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
22.
Anderson
BJ
,
Holford
NH
:
Mechanistic basis of using body size and maturation to predict clearance in humans.
Drug Metab Pharmacokinet
2009
;
24
:
25
36
23.
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
24.
Glass
PS
,
Shafer
S
,
Reves
JG
:
Intravenous drug delivery systems, Miller’s Anesthesia
. Edited by
Miller
RD
.
Philadelphia, Pennsylvania, Elsevier
,
Churchill Livingstone
,
2005
, pp
439
480
25.
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
26.
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
27.
Anderson
BJ
,
Holford
NH
,
Woollard
GA
,
Chan
PL
:
Paracetamol plasma and cerebrospinal fluid pharmacokinetics in children.
Br J Clin Pharmacol
1998
;
46
:
237
43
28.
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