Disposition of Remifentanil in Obesity: A New Pharmacokinetic Model Incorporating the Influence of Body Mass

REMIFENTANIL is an esterase metabolized µ opioid agonist in widespread clinical use internationally, often as part of total intravenous anesthetic techniques. The clinical pharmacology of remifentanil has been exhaustively investigated, perhaps more so than any previously introduced fentanyl congener. Population mixed-effects models describing the pharmacokinetic (PK) and pharmacodynamic (PD) behaviors of remifentanil in quantitative terms have been developed and validated.^1–4  The influences of numerous covariate effects, such as age, sex, kidney function, and hepatic function, have also been characterized (some of these covariate effects are incorporated into PK–PD models that are widely applied clinically).^3,5,6  The synergistic interaction of remifentanil when combined with propofol or inhaled anesthetics has also been modeled and evaluated in the clinical domain.^7–11 

Despite this substantial collection of PK–PD models characterizing the disposition, effects, and interactions of remifentanil, a clear unmet need in our understanding of the clinical pharmacology of remifentanil relates to the influence of obesity on remifentanil PKs. Previous investigations addressing this issue suffer from some important limitations despite using rigorous data-gathering methods and advanced model-building techniques. The study by Minto et al.,³  the most widely used model in target-controlled infusion (TCI) systems, was methodologically flawed as it relates to the obesity covariate, because the equations (i.e., the James’ Equations) used to calculate lean body mass in this study are now known to be seriously flawed, particularly in very obese patients.¹²  More importantly, the Minto model was estimated from volunteer subjects that did not include obese subjects. The study by Egan et al.,¹³  although focused entirely on the influence of obesity on remifentanil PKs, must be regarded as preliminary and incomplete, because the study was relatively small and, although the study compared obese and lean groups, the study population did not include large numbers of morbidly obese patients. Studies by La Colla et al.^14,15  also address remifentanil PKs in the obese patient population, but these studies offer workarounds to address the flawed lean body mass calculations of the study by Minto et al.³  rather than a new PK model that clearly incorporates body weight indices into the parameters of a PK model. The work by La Colla et al.¹⁵  clearly demonstrates that the performance of the PK parameters from Minto et al.³  lead to a large bias when applied to morbidly obese patients.

The relevance of the current study is in part a function of the international obesity epidemic. Given the high prevalence of obesity around the world and the high incidence of obesity-associated disease, obese patients frequently present for anesthesia and operations, sometimes to treat their obesity by a bariatric procedure.^16–19  Understanding the influence of obesity on the disposition of anesthetic drugs is therefore an important question in contemporary anesthesia practice.²⁰  The aging of the worldwide population in many developed countries and the high incidence of age-related comorbidities make understanding the influence of age on anesthetic drug disposition and effects similarly important, although the age issue has already been explored in detail, especially for advanced age.

The aim of the current investigation was to develop a new population PK model for remifentanil that would adequately characterize the influence of body weight (among other covariates, e.g., age) on the disposition of remifentanil in the general adult population. We hypothesized that age and various indices of body mass would be important covariates in the new model. Our practical aim was to build an updated PK model for remifentanil suitable for general adult use, including obese patients.

The data set for model building was derived from nine previously published studies by Egan et al.,^1,2,13,21  Minto et al.,³  Drover and Lemmens,²²  Mertens et al.,⁴  Kern et al.,²³  and La Colla et al.¹⁵  All of the studies were approved by local human-subject institutional review boards, and all of the subjects gave informed consent (see individual studies for details). In total, there were 4,455 remifentanil concentration measurements from 229 subjects (107 patients and 122 volunteers); all of the subjects and data points from the nine previously published studies were included in the data set for analysis. The studies included remifentanil administration by bolus injection, continuous infusion, and TCI. Blood sampling ranged from 2 to 523 min. The measured plasma concentrations ranged from 0.054 to 245.4 ng/ml. Details about each of the nine studies composing the data set are presented in tables 1 and 2.

Population PK analyses were performed using NONMEM (version 7.2, Icon Development Solutions, USA). Three-compartment models with linear PKs were fitted using ADVAN6 subroutines and the first-order conditional estimation with interaction procedure.

A log normal variance model was used to describe the interindividual variability of the remifentanil PK parameters according to the following equation:

where θ_TV is a population mean value for the PK parameters of the population, θ_i is the individual post hoc estimate for the PK parameter in the ith subject, and η_i is a random variable that represents the difference between individual (θ_i) and population (θ_TV) values that is a normally distributed random variable with a mean of 0 and a variance ω².

For the intraindividual variability that describes the residual errors, combined error models were tested using the following equation:

where C_ij is the j th plasma concentration measured in the i th subject; C_pred,ij is the j th concentration predicted based on the model parameters, dosing regimen, and sampling time in the i th patient; and ε_pro,ij and ε_add,ij are the coefficient of proportional residual error and coefficient of additive residual error terms, respectively, which are normally distributed random variables with a mean of 0 and variance of σ².

The potential covariates affecting model structural parameters were explored for body weight, age, sex, body mass index (BMI), and fat-free mass (FFM) with the formula by Janmahasatian et al.²⁴  FFM was calculated from sex, weight (in kilograms), and BMI (in kilograms per meter squared):

To explore the covariate-parameter relationships that could further explain interindividual variability, the estimated parameters obtained from a three-compartment model without any covariates were plotted independently against body weight, age, sex, BMI, and FFM. Those graphical relationships were examined visually to identify potentially important covariate effects. Promising covariates were then entered into the base model sequentially (with either a linear or power function centered on the median of the covariate depending on the exploratory plots).

The process for including covariate effects within the mixed-effects structural model proceeded as advocated by Jonsson and Karlsson.²⁵  Starting from the base model, potentially important covariate-parameter relationships were individually tested in a stepwise fashion building the model forward (i.e., the forward inclusion process). The most influential covariates in terms of the NONMEM objective function value (OFV) reduction were iteratively retained in the model. The minimal OFV reduction (equal to minus twice the log likelihood) statistically justifying inclusion of a candidate covariate effect into the model was determined by a P value of less than 0.01, that is, a decrease in the OFV of at least 6.63 points for one additional parameter (using a chi-square distribution, 1 degree of freedom).

When no more covariate effects were statistically justified for inclusion into the model, the so-called final model was defined. Thereafter, a backward elimination procedure was performed in which each covariate was removed in turn from the final model, and the difference in OFV between the final and each reduced model was examined. An increase in OFV of 10.83 (P < 0.001) was required to retain the covariate in the final model (i.e., the threshold for the change in OFV during the backward elimination procedure was set more parsimoniously than during the forward inclusion of covariates).

Goodness of fit for all models was visually inspected with plots of observed versus predicted values, looking for visual evidence of inaccuracy and bias. The conditional weighted residuals (CWRES) of the PK models were plotted as a function of time to assess model misspecification.²⁶  Predictive performance and log likelihood profiles were examined, facilitated by Wings for NONMEM, xpose4 (version 4.0), and fit4NM run on the R statistical software package (version 2.13.1, the R Foundation for Statistical Computing, Austria).²⁷  As an internal validation procedure, a nonparametric, bootstrap analysis was conducted.

The predictive performance of the models was evaluated numerically by computing the performance error (PE).²⁸  For each blood sample the PE was calculated as follows:

where C_p,ij is the predicted blood remifentanil concentration in j th sample from the i th patient, and C_m,ij is the measured blood concentration of remifentanil in that sample. Subsequently, the median prediction error (MDPE) and median absolute prediction error (MDAPE), which indicate bias and inaccuracy, respectively, were calculated in a pooled data method as follows:

where MDPE_i and MDAPE_i are median PE and median absolute PE of i th individual, PE_ij is the PE in j th sample from i th individual, N_i is the number of observation in i th individual, and M is the number of individuals.

The robustness of the MDPE and MDAPE calculations were evaluated with a 2-fold cross-validation procedure. With this method, the data set was randomly divided into two parts: the first arm was used to estimate parameters and the second arm was used to predict concentrations with the parameters derived from the first arm. This process was repeated after exchanging the arms. The prediction results of the two arms are then combined and treated as one data set. This process was repeated 10 times with different random partitioning of the data set. All of the results were pooled, and MDPE and MDAPE were calculated.

To evaluate the robustness of the final model, a nonparametric bootstrap analysis was performed as an internal validation procedure.²⁹  Two-thousand bootstrap data set replicates were generated randomly by resampling with replacement (i.e., each data set included 229 subjects selected randomly; in each data set, some subjects might be randomly represented more than once and some might be randomly omitted). Model parameters for each of these data sets were estimated using NONMEM. Median parameter values and the 2.5 to 97.5 percentiles from the bootstrap procedure were compared with those of the final model parameter estimates. If the final model parameter estimates were not significantly different from the results of the bootstrap procedure, the model was considered stable.

Log-likelihood profiles (LLPs) for each of the final model parameters were computed as a means of evaluating uncertainty in the parameters. LLPs are computed by fixing a given model parameter at values around the final model estimate and then rerunning the model to obtain new OFVs. Plots of the fixed parameter values versus the OFVs (i.e., the LLPs) are intended to reveal whether there is a problem with parameter identification. The LLP determines the parameter values on either side of the final estimate that produce a deterioration in the OFVs of 3.84. This relies on the assumption that the difference of –2LL follows a chi-square distribution and that when this value changes by 3.84, the parameter value is at the 95% confidence limit.

Deterministic computer simulations using the covariate adjusted population PK model were performed to illustrate its clinical implications. First, we explored the time course of remifentanil plasma concentration (Cp) when the same absolute dose (i.e., not weight adjusted) is administered intravenously to 25-yr-old subjects with different body weights (i.e., 75 and 150 kg) and a height of 175 cm, and also to 75-kg subjects with different ages (i.e., 25 and 75 yr old) and a height of 175 cm. The simulated dosing regimens were a 50-µg bolus injection and a 15-µg/min continuous infusion for 60 min.

Next, to illustrate the difference in remifentanil dosing requirements for various body weights (ranging from 50 to 200 kg in a 50-yr-old subject, height of 175 cm) and ages (ranging from 25 to 75 yr old in a 75-kg subject, height 175 cm), we performed TCI simulations, which targeted a Cp of 5 ng/ml for 30 min. We also simulated the cumulative dosage required to achieve a Cp of 5 ng/ml for various ages and body weights for up to 300 min.

Finally, because the final covariate added to the model (relating FFM to volume of distribution of the second compartment [V2]) was associated with the smallest change in the NONMEM OFV, we explored through simulation the influence of differing FFM values (i.e., 45, 55, and 65) on the predicted plasma concentrations (and doses) using the final model (wherein FFM is a covariate related to V2) and the penultimate model (i.e., model 8 in table 3, where FFM is not a covariate for any parameter). We simulated bolus, infusion, and TCI administration to a 50-yr-old, 75-kg patient (with arbitrarily different FFM values for the final model, FFM of 45, 55, and 65). The simulated remifentanil bolus was 50 µg. The simulated remifentanil infusion was 20 µg/min for 20 min. The simulated TCI regimen was a 30-min infusion targeted to 5 ng/ml.

All of the simulations were implemented in PKPD Tools (Minto and Schnider, http://www.pkpdtools.com, accessed March 21, 2016).

The study population demographics for age and body weight are presented graphically in figure 1 (see also table 1). The demographic characteristics confirm a heterogeneous population, particularly regarding the covariates of interest (i.e., body weight and age); BMI ranges from 16.1 to 73.7 and age from 20 to 85 yr. The population includes patients with extremely high body weight and BMI from the subjects of La Colla et al.¹⁵  and Egan et al.¹³  Numerous older subjects are also represented, particularly in the Minto et al.³  and Egan et al.²¹  data sets.

Not only are the demographic characteristics of the study population heterogeneous, but the dosing schemes are also highly variable in terms of method of administration and total dose. The dosing scheme, concentration range, sampling site, number of samples, sampling period, and volunteer/patient distribution for each study are shown in table 2. The dosing schemes included high-dose bolus injection, fixed-rate continuous infusion, and TCI infusion. All of the concentration measurements were made from arterial blood. The raw PK data, categorized by method of administration, are presented in figure 2.

The mixed-effects modeling analyses showed that the PKs were best described by a three-compartment model (clearance [CL1], intercompartmental clearances [CL2 and CL3], volume of the central compartment [V1], and volumes of peripheral compartments [V2 and V3]). The process of model development is presented in summary form in table 3, with OFVs and numerical performance measures for a two-compartment model, the base model, and for selected covariate-adjusted candidate models (hundreds of candidate models were tested). The parameter estimates and interindividual variability for the base model and the covariate adjusted final model are presented in table 4 (bootstrap median parameter values with CIs shown for the final model).

In the process of covariate model building, age, weight, and FFM were discovered to have an important influence on several model parameters. Sex had no influence on any model parameter (once body weight was incorporated into the model). Body weight was related to clearance (CL1) and central volume (V1) as a power function with a positive exponent smaller than one. Central clearance (CL1), fast intercompartmental clearance (CL2), and rapidly equilibrating and slowly equilibrating compartment volumes (V2 and V3) exhibited linear correlations with age. FFM was modeled into the fast compartment volume (V2) as a power function with a positive exponent smaller than one. The relationships between parameters and their related covariates are presented as Supplemental Digital Content 1 (https://links.lww.com/ALN/B427).

The observed versus predicted remifentanil concentrations for the final population model and the post hoc individual models are presented in figure 3. Similar plots for the penultimate model (i.e., model 8 in table 3) are presented as Supplemental Digital Content 2 (https://links.lww.com/ALN/B428). The CWRES as a function of time for the final population PK model are presented in figure 4. The percentage of CWRES outside of the ± 2 range is 4.93%. The MDPE and MDAPE with CIs are presented in table 5 (and also for the Minto and La Colla models). Compared with the Minto and La Colla models, the final model has minimal bias (for MDPE, the Minto and La Colla model CI excludes 0, an indicator of bias). The results of the 2-fold cross-validation procedure are also presented in table 5; the values differ minimally with the original values. Plots of the post hoc ηs for CL1 versus weight and CL1 versus age for the base and final models are presented as Supplemental Digital Content 3 (https://links.lww.com/ALN/B429).

The coefficients of variation and the CIs (according to the bootstrap procedure) of the parameter estimates for the final model are presented in table 4. The median values produced from the bootstrap procedure are close to the estimated parameters. Reduction of the coefficients of variation of parameters was achieved in the final model compared with the base model except for the slow intercompartmental clearance.

The LLP plots are presented in figure 5 for all of the estimated parameters in the final model. This LLP analysis suggests that the parameters were estimated with adequate precision (i.e., no major problems with parameter identification).

The same-dose simulations are presented in figure 6. Given the same dose, an obese patient is expected to achieve modestly lower remifentanil plasma concentrations than a lean patient of the same age (fig. 6A). Conversely, an older patient is expected to achieve substantially higher plasma concentrations than a younger patient of the same body weight after an identical dose (fig. 6B).

The changes in dosing requirements over time for a range of body weights and ages when targeting a remifentanil Cp of 5 ng/ml by TCI are presented in figure 7 (fig. 7A for weight and 7B for age). Because of the relatively high clearance of remifentanil, the infusion rate required stabilizes approximately 25 min after beginning the infusions. The increased dose requirements for the younger and heavier patients are evident.

The cumulative dose requirement to maintain a Cp of 5 ng/ml by TCI simulation is presented in figure 8. For a 5-h infusion, compared with the younger-leaner patient (i.e., 25 yr old, 75 kg), the heavier (150 kg) 25-yr-old’s cumulative dosage requirement is approximately 125% of the younger-leaner patient’s dose. Conversely, the older (75 yr old) 75-kg patient’s cumulative dosage requirement is approximately 75% of the younger-leaner patient’s dose. Table 6 presents how dosing requirements vary by body weight for infusions of up to 2 h (for a 175-cm 25 yr old).

The simulations comparing the penultimate model (model 8 in table 3) with the final model for differing FFM values in an otherwise identical patient are presented as Supplemental Digital Content 4 (https://links.lww.com/ALN/B430). As expected, the differences between the model predictions are slight. These simulations must be interpreted with caution, because assuming such a wide range of FFM in a 75-kg patient is somewhat artificial.

The aim of the current investigation was to develop a new population PK model for remifentanil that would adequately characterize the influence of body weight (among other covariates, e.g., age) on the disposition of remifentanil in the general adult population. We hypothesized that age and various indices of body mass would be important covariates in the new model. This hypothesis was confirmed.

That remifentanil PKs are influenced by obesity and age is of course not new. This work extends the findings of previous investigators, such as Minto et al.,³  La Colla et al.,¹⁵  and Egan et al.,^13,21  by including data from a larger number of subjects with a broader range of body weight and age and a more diverse array of drug infusion schemes, including TCIs. Compared with the data sets for typical descriptive PK studies in anesthesiology, a model derived from more than 225 subjects and nearly 4,500 arterial blood samples is a large, comprehensive study. Overcoming the limitations associated with the faulty James’ equations for the computation of lean body mass constitutes a simple but important advance of the existing knowledge.

From a practical perspective, PK models have at least two important applications.³⁰  First, the models can be the basis for PK simulations to address clinically important questions regarding the optimal dosage scheme for a given patient (i.e., posologic optimization). An extension of this application is the incorporation of the models into clinical pharmacology display systems for clinical and educational use.^31,32  Second, PK models can be incorporated into TCI pumps for the automated administration of intravenous anesthetics, enabling practice in the concentration domain.³³ 

In terms of the practical treatment of obese patients, the application of this new PK model is straightforward. Using total body weight (TBW) as a weight-proportional approach to calculate the drug administration scheme in severely obese patients is suboptimal. Compared with lean subjects, severely obese patients do indeed require more remifentanil to achieve a specified remifentanil target concentration. However, in considerably obese patients, the dosage increase does not correlate linearly with their TBW (in fact, it is not even a close approximation of a linear increase). This is consistent with existing information about the influence of body weight on remifentanil disposition.^13,15  According to the current study, although obesity does indeed influence remifentanil PKs, it is not necessary to relate basic or direct obesity indices, such as BMI, as a covariate to any of the PK model parameters. A model that incorporates TBW, fat-free mass, and age as covariates is sufficiently accurate, even for obese patients.

From a theoretical perspective, the scientific underpinnings of this relationship between PK model parameters and body weight are a function of well-characterized anatomic realities. In obese patients, as body mass increases, fatty tissue increases more than lean body mass; that is, there is a nonlinear relationship between TBW and lean body mass.³⁴  Lean body mass, especially as it relates to clearance and metabolism, is more relevant pharmacologically than fatty tissue, although fat is thought to be an important reservoir for drug distribution for fat-soluble drugs.³⁵  Because remifentanil is metabolized in blood and tissue by nonspecific esterases that are presumably widely expressed in the body,³⁶  it is conceivable that TBW is perhaps more relevant to the metabolism of remifentanil than drugs that require delivery to a metabolic organ like the liver; that TBW is a covariate in the final model is consistent with this speculative assertion.

Recent advances in understanding the disposition of propofol in the obese patient population are instructive to put the current findings in proper context. Cortínez et al.³⁷  reported an allometric scaling relationship between propofol elimination (and intercompartmental clearances) and TBW (i.e., TBW^3/4). In a much larger and comprehensive study of numerous propofol data sets pooled together, Eleveld et al.³⁰  confirmed the use of allometric scaling of propofol clearance parameters. In contrast, in a clinical study, Ingrande et al.³⁸  concluded that lean body weight was the optimal basis for the calculation of the propofol induction dose. A sophisticated simulation study by McLeay et al.³⁹  also concluded that lean body weight was a useful descriptor for the determination of propofol dosing schemes. In general, although these studies have nuanced differences regarding how best to include the influence of obesity into an optimized PK model for propofol, they all confirm the fundamental principle that TBW is not an optimal scalar for the calculation of propofol dosing schemes. The current study findings about remifentanil are in general agreement with these observations for propofol.

Investigators in this arena have focused considerable attention on allometric scaling of PK parameters to body weight. In general terms, allometry is the study of the relationship between body size and body anatomical and physiologic properties. The theoretical foundation of allometric scaling in PKs is loosely based on Kleiber’s law, which relates the body size of numerous species to metabolic rate using a ¾ power function (i.e., metabolic rate is proportional to TBW^3/4). Extending this idea into PKs, numerous investigators have advocated scaling PK clearance parameters to TBW^3/4.^12,40  As noted, this approach has been successful for propofol.^30,37  There is controversy in the literature about how best to implement some form of allometric scaling in PK models.⁴¹  The current study took a hybrid approach, estimating the exponent of the power scaling as part of the model building (i.e., rather than using a ¾ power function). In other words, the current study relied exclusively on the data to determine the structure of the model, including the linear or power scaling of the covariates, rather than applying these allometric scaling principles from the outset. This approach is consistent with a recent editorial cautioning blind adherence to allometric PK model building strategies when the raw data do not support it.⁴² 

Several limitations of this investigation deserve emphasis. The most obvious limitation is that the study does not include a PD component. PKs cannot be considered in isolation, so the study conclusions must be understood in the context of existing knowledge regarding remifentanil PDs. Remifentanil PDs are not thought to be influenced by obesity, and the influence of advancing age on remifentanil PDs is well characterized.³  However, how to predict effect-site concentrations using the PK model of the current study is a conundrum, because we do not have a ke0 parameter (elimination rate constant from the effect-site) estimated from this data set. Using the time of peak effect after bolus administration is one approach to this problem.⁴³ 

Another significant limitation of this study relates to the remifentanil blood samples that were below the quantitation limit (BQL) of the assay. It is well known that some methods used to handle the BQL values can introduce bias into model parameter estimates, especially where the data sets are sparsely sampled.^44,45  Because the current study relied on data from numerous investigators and studies, some of which were published many years ago, it was not possible for us to determine with confidence how many samples were BQL, although we believe that, relative to the total number of samples, the number of BQL samples was very small (perhaps less than 2%) and that all of the studies had a similar limit of quantitation. We elected to use all of the data from the published studies because we did not have a rational basis to censor the data further. There is no question that combining several studies that used different analytical chemistry techniques and laboratories may have influenced the final model because of assay measurement error and that this issue is especially relevant to BQL samples (or samples that are at or near the limit of quantitation).

Other limitations include the absence of pediatric data, making the model unsuitable for use in younger patients. Other covariate effects that might be of interest were also not examined, including subject type (volunteer vs. patient). Effects for these covariates have been described for propofol.³⁰  Although in preliminary models we identified a very modest covariate effect for subject type (volunteer vs. patient), we did not pursue this covariate effect in the final model in part because all of the obese subjects were patients, and we were concerned about this confounding influence (our main goal was to focus on the influence of obesity).

The next steps in this line of investigation obviously include prospective validation of the new PK model. Assessment of the predictive performance in various subgroups would also be enlightening. Finally, there is considerable work to do in illustrating the clinical implications of the new PK model through simulation.

The general question of how to handle dosing calculations in very obese patients has important implications for drug development in all therapeutic areas. Given the prevalence of obesity in the modern patient population, particularly in developed countries like the United States, the development of a drug label should include information about how dosing schemes should be modified with increasing body weight.^46,47  This is perhaps even more important in anesthesiology, where drugs of low therapeutic index are common. The current study provides a rational basis to address these obesity-related posologic issues for remifentanil.

In summary, we report a new PK model for remifentanil estimated from a large, diverse data set that adequately characterizes the influence of body weight and age on the disposition of remifentanil in a general adult population. The model is suitable for use in TCI systems and pharmacologic simulation software. The model requires prospective validation before widespread use.

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

The authors declare no competing interests.

The Remifentanil Pharmacokinetics in Obesity Investigators for this study are as follows: Charles F. Minto, M.D., Ph.D., Luca La Colla, M.D., David R. Drover, M.D., Jaap Vuyk, M.D. Ph.D., Martijn Mertens, M.D.