An Allometric Model of Remifentanil Pharmacokinetics and Pharmacodynamics

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 displays¹  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.²  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 children³  and adults.⁴  Combining these data might lead to a more general applicable PK-PD model for remifentanil. A similar approach was taken by Eleveld et al.,⁵  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,⁹  more recently it has been systematically applied to PK models,¹⁰  and it has been used for remifentanil.¹¹  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.

We obtained remifentanil PK and PD data from three studies: the aforementioned study in adults by Minto et al.,²  a study in children by Ross et al.,³  and a study in adults with TCI by Mertens et al.⁴  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.

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,¹²  version 2.14.1. Simulations of TCI using the models were performed using PK-PD Tools for Excel^© (www.pkpdtools.com).

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.¹⁰  We also considered the fat-free mass (FFM) predictor described by Al-Sallami et al.¹³  as a body-size descriptor that uses a maturation model to extend the Janmahasatian et al.¹⁴  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,¹⁷  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 ω² is the estimated parameter population variance.

To quantify the PK predictive performance for an observation, we calculated the performance error¹⁸  (PE_PK) and absolute performance error (APE_PK) as follows.

For PD measures, only the Minto study data provided PD observations that were the spectral-edge frequency (SEF) of the electroencephalograph.¹⁹  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 PE_PK (MdPE_PK), median APE_PK (MdAPE_PK), median PE_PD (MdPE_PD), and median APE_PD (MdAPE_PD). 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 MdAPE_PK and MdAPE_PD across these subgroups. The subgroups were chosen to approximately parallel PK models in the literature.

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,²  the Egan model,²⁰  and the La Colla model.²¹  We also considered the Rigby-Jones model,¹¹  which is an allometric PK model developed from data from children.

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.

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.

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 allometry⁵  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/V2_ref)^0.75 and (V3/V3_ref)^0.75, respectively instead of the usual allometric approach of scaling with TBW, i.e., with (TBW/TBW_ref)^0.75.

We tested an exponential age covariate relationship centered at 35 yr¹⁷  for all parameters (Model 5). This involves multiplying each parameter by a function of the form f_aging(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 f_sigmoid(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 MdAPE_PK 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, https://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 (https://links.lww.com/ALN/B422). The likelihood profiles (Supplemental Digital Content 3, https://links.lww.com/ALN/B423) and the relationships between η and the covariates (Supplemental Digital Content 4, https://links.lww.com/ALN/B424) and between η1 to η6 (Supplemental Digital Content 5, https://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 V1_ref, V2_ref, V3_ref, CL_ref, Q2_ref, and Q3_ref are the estimated compartmental volumes and clearances for the reference individual, Θ₁ defines the weight at which 50% of maturation of CL is complete, Θ₂ defines the change in V1, Q2, and Q3 with age, Θ₃ defines the change in V2 and CL with age, Θ₄ defines the change in V3 with age, Θ₅ defines the change in V2, CL, and Q2 for females aged 12 to 45 yr, and Θ₆ 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.

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.

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.

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.

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 (https://links.lww.com/ALN/B422). The likelihood profiles can be found in Supplemental Digital Content 6 (https://links.lww.com/ALN/B426). The summarized equations of the final models are as follows.

Symbols E0_ref, Emax_ref, Ce50_ref, γ_ref, ke0_ref, and Θ₁ 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.

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.

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.

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.²⁵  Additional data are needed to determine the performance of our model in the obese.

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.

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 elimination²⁶  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-life²⁷ ). 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.

As reported for propofol,⁵  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 Holford¹⁰  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.

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.

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

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.

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.

Support was provided solely from institutional and/or departmental sources.

The authors declare no competing interests.