A Model of the Ventilatory Depressant Potency of Remifentanil in the Non–steady State

The study was approved by the Stanford University Institutional Review Board (Stanford, California). Written informed consent was obtained from each subject. The reported data are a subset from a study of propofol and remifentanil aimed at identification of either drug's ventilatory depressant effects, the pharmacokinetic interaction between both drugs, and the interaction of both drugs with regard to suppression of quantal responses to central nervous system stimulation and electroencephalographic effects.

We studied five male and five female healthy volunteers (age, 33.5 [23–43] yr, weight, 69.3 [50–100] kg). All volunteers received a physical examination, laboratory tests (complete blood cell count, blood chemistries), and an electrocardiogram.

All volunteers were studied after fasting for at least 6 h. After arrival in the operating room, standard monitoring (noninvasive blood pressure, electrocardiography, and pulse oximetry) was established; one arterial cannula (radial artery of the nondominant hand) and two intravenous cannulae (both forearms) were inserted. The volunteers were supplied with a continuous positive airway pressure mask connected to a pressure differential spirometer/sidestream gas analyzer (D-lite flow sensor/gas sampler, AS/3; Datex-Ohmeda, Helsinki, Finland). Bispectral Index was recorded with an Aspect 1000 electroencephalographic monitor (BIS® monitor Aspect Medical, Natick, MA). Before the study, the pressure differential spirometer underwent three-point calibration (500, 1,000, 1,500 ml) with a 3-l calibration syringe (Hans Rudolph, Kansas City, MO). The gas analyzer underwent two-point calibration with gas mixtures containing 4% and 8% CO². Heart rate, blood pressure, ventilatory rate, tidal volume, MV, and inspiratory/expiratory oxygen and carbon dioxide were recorded every 5 s using the software program Collect (Datex, Helsinki, Finland).

After a 5-min resting and equilibration period, an arterial blood sample for determination of arterial carbon dioxide tension (Paco²) was drawn.

Two anesthesia ventilation bags (volume, 2.3 l each) connected with a Y piece were filled with oxygen and connected to the D-lite ventilatory sensor. The volunteers breathed from this reservoir and therefore rebreathed their exhaled carbon dioxide. Contrary to the classic Read design, 8the bags contained no carbon dioxide when rebreathing was initiated. At the volunteers’ request, the rebreathing bags, but not the D-lite sensor, were removed (on average after 5 min), enabling us to obtain blood gas and volume measurements during recovery to normal resting ventilation. At the start of and during and after the rebreathing part of the study, arterial blood samples were drawn every 1–2 min for the determination of Paco². After stabilization of MV at baseline levels, this part of the study was terminated.

After determining individual carbon dioxide responses, remifentanil was administered via  target-controlled infusion with a Harvard infusion pump (Harvard Clinical Technology, Inc., South Natick, MA) driven by STANPUMP #running on a commercially available laptop computer. The remifentanil pharmacokinetic parameters were the covariate adjusted set reported by Minto et al.  9Remifentanil was administered in ascending steps targeting plasma concentrations until end-tidal carbon dioxide (Peco²) exceeded 65 mmHg and/or apnea periods of more than 60 s occurred. Thereafter, the remifentanil concentration was allowed to passively decrease to 1 ng/ml. Table 1displays the target concentrations steps for each volunteer. During and after the infusion, frequent arterial blood samples for determination of remifentanil concentrations and Paco²were drawn. Figure 1displays a typical example of this administration schedule and the corresponding Paco²values.

Blood sampling was timed based on pharmacokinetic, pharmacodynamic, and efficiency considerations (high-resolution sampling during periods of rapidly changing concentrations). A blank sample was drawn after insertion of the arterial cannula. Blood samples were drawn 2, 5, 10, and 15 min from the start of the infusion. For every further step up, one sample was drawn immediately before changing the target concentration. During the passive decrease down to 1 ng/ml for remifentanil, samples were drawn at 2, 4, 6, 8, 10, 15, and 30 min after the last change in the target concentration. Citric acid was added to the samples immediately after blood was drawn. Samples were stored at −20°C until assaying (for assay method see Bouillon et al.  10).

All arterial blood samples drawn for analysis of Paco²were stored on ice immediately after the drawing of blood and were analyzed within 20 min with a portable blood gas analyzer (i-stat Corporation, East Windsor, NJ). All volume measurements were converted to standard temperature and pressure before entering further calculations. Because Paco²measurements were at best available every minute, MV measurements every 5 s and the further analysis required data pairs of Pco²and MV matched by time, Peco²values were used for calculation of the carbon dioxide response. To check for the validity of this approach, the average difference between Peco²and Paco²and the corresponding 95% confidence interval was calculated. Paco²measurements were used for calculation of drug-induced ventilatory depression.

The pharmacokinetic model has been described previously. 10In brief, the arterial concentration time course of remifentanil was adequately predicted by a two-compartment model with the following parameters (typical value ± coefficient of variation), V¹= 13.50 l, V²= 8.64 l ± 27.4%, Cl¹= 2.57 l/min ± 25.3%, Cl²= 1.31 l/min ± 7.7%. The individual post hoc  Bayesian predictions of this model were used for the pharmacodynamic calculations.

Proportional and exponential models were used to describe the interindividual variability of the pharmacodynamic parameters:

where θ⁽ⁱ⁾refers to the individual value of the respective parameter in the i  th individual, θ^(TV)is the typical value of the respective parameter in the population, and η varies randomly between individuals with mean zero and variance of Ω. 2 

An additive (constant SD) error model was chosen for modeling residual variability of both MV and Paco²:

DV^obsrefers to the observed value of the dependent variable (MV, V̇^alv); DV^exprefers to the value predicted based on dose, time, and the individual pharmacokinetic and pharmacodynamic parameters. ε is a normally distributed random variable with mean zero and variance σ. 2 

The program system NONMEM, version V with the “First Order Conditional Estimation” method and “η-ε interaction,” was used for all model fits and empirical Bayesian estimation of the individual parameters. 11 

Decisions between different models were made using the log likelihood test (P < 0.01). Model misspecification was checked for by plotting the predicted against the measured values of the dependent variable.

Minute ventilation, MV, was modeled as a function of end-tidal carbon dioxide, Peco², the independent variable. Hysteresis in the MV versus  Peco²relation was modeled using an effect compartment for carbon dioxide:

where Pecco²is the partial pressure of carbon dioxide at the effect site (“biophase”; mmHg) and k^e0,CO2is the first-order equilibration constant between arterial and effect-site Pco².

Before rebreathing, the system is at steady state, and baseline carbon dioxide, Pecco²(0), equals Peco²(0).

Although the relation between MV and Pecco²above the metabolic hyperbola (the steady state relation between alveolar ventilation and Paco²for a given carbon dioxide production and inspiratory fraction of carbon dioxide) can be well described by a straight line, 12,13the shape of the curve changes near the hyperbola. We speculated that this change persists, and is even exaggerated, if ventilation is depressed below baseline. To account for our hypothesized shape of this relationship, Pecco²was used as an independent variable of a nonlinear expression 6:

where MV(Pecco²) is the MV depending on Pecco²(l/min); MV(0) is the baseline MV (l/min); Pecco²(0,t) is the partial pressure of carbon dioxide at the effect site, baseline, and time t (mmHg); and F is the gain determining the change of MV for a given ratio of Pecco²(t) and Pecco²(0). For reasons of completeness, a linear carbon dioxide response was also calculated as:

where SL is the slope of the carbon dioxide response curve. For each individual, baseline Pecco²was fixed to the measured value of baseline Peco².

Because changes in Paco²during drug-induced ventilatory depression are not as fast as those during carbon dioxide rebreathing and Paco²is less artifact sensitive than Peco²(mask fit, shallow breathing), Paco²was used as the dependent variable for modeling of remifentanil-induced ventilatory depression. Although described previously, 6the essential steps of the modeling approach are repeated here.

Changes of partial pressures of a gas in the body over time can be computed with mass balance equations. For a one-compartment model with constant input (carbon dioxide production) and constant output (carbon dioxide elimination) under baseline steady state conditions, the change of amount of carbon dioxide over time can be expressed as:

where V^dco²is the apparent volume of distribution of carbon dioxide (l), Paco²is the arterial partial pressure of carbon dioxide (mmHg), kⁱⁿis the production rate of carbon dioxide (l/min), k^outis the elimination rate of carbon dioxide (l/min), and 760 is the atmospheric pressure at sea level (mmHg). k^out(t) can also be expressed as the product of alveolar ventilation (l/min) and the current Paco²divided by the barometric pressure, yielding:

Under the assumption that the production rate of carbon dioxide is always equal to the baseline elimination rate, the production rate can be substituted by the product of the baseline value of the normalized Paco²and alveolar ventilation and becomes a constant:

Rearranging for the change of Paco², the dependent variable, over time yields:

Introduction of a hypothetical effect compartment for remifentanil was attempted according to the following equation but did not improve the fit of the model to the data:

where Cp is the drug concentration in plasma calculated from the individual dosing histories and pharmacokinetic parameters, Ce is the drug concentration in the effect compartment, and k^e0is the first-order rate constant governing the transfer of drug out of the effect compartment.

The combined inhibitory effect of remifentanil (plasma concentration) and the stimulatory effect of Pecco²on alveolar ventilation was then expressed as product of a fractional sigmoid Emax model and the nonlinear term for carbon dioxide response: where V̇^alv(0) refers to baseline alveolar ventilation, c^prefers to the plasma concentration of remifentanil, and C⁵⁰refers to the remifentanil concentration at which V̇^alvand therefore k^outwill be decreased to 50% of the value in the absence of remifentanil, for unchanged Pecco². F was calculated from the individual carbon dioxide response curves. This equation also yields alveolar ventilation normalized to baseline (divide both sides by V̇^alv(0)) and can therefore be used for predictions of the time course of ventilation after the administration of drugs, even in the absence of measured ventilatory data.

After insertion into the mass balance equation, the final equation to describe opioid induced hypercapnia can be obtained:

As an alternative approach, we used the power function advanced by Dahan et al.  7:

One volunteer was excluded from the analysis because of a very low baseline Paco², presumably caused by anxiety; in one volunteer, carbon dioxide rebreathing failed as a result of equipment problems. No volunteer was more than lightly sedated, as judged from both clinical observation and electroencephalographic data.

The parameters of a nonlinear model characterizing the influence of carbon dioxide on MV including effect compartment equilibration are summarized in table 2. Table 3shows the respective values for a linear model. Applied to the entire range of measurements, a nonlinear model described the data better than a single linear model, with a difference in objective function of 388 (P < 0.001). Therefore the predictions of the linear model were not plotted, and its parameters were not used for further calculations.

The population and Bayesian predictions and goodness of fit of the nonlinear model are summarized in figure 2. The upper panel shows the carbon dioxide response curves of the volunteers (Bayesian predictions of MV vs.  Pco²in the effect compartment [Pecco²]). The lower panel shows the Bayesian predictions of MV versus  the measured MV. Because the data points are symmetrically distributed around the line of identity, the effect compartment model and nonlinear model for carbon dioxide response adequately capture the relation of MV and carbon dioxide in the range of measurements. The Bayesian predictions of F, the gain of the nonlinear carbon dioxide response curves, and k^e0,CO2were used for the calculation of remifentanil-induced ventilatory depression.

The remaining parameters characterizing the indirect response model describing remifentanil-induced ventilatory depression are summarized in table 4. Note that k^el,CO2is the elimination constant of carbon dioxide and completely independent of drug action. Paco²(0) was measured. The sigmoid Emax pharmacodynamic model (equation 12) resulted in a 25-point improvement in log likelihood over the power pharmacodynamic model (equation 13), which was therefore not considered further.

Figure 3displays the population and Bayesian prediction of Paco²for one volunteer (top ) and the goodness of fit of both the population and Bayesian estimates for all volunteers (bottom ). Although the model is not misspecified, as can be seen from the good fit of the individual post hoc  Bayesian predictions, the amount of scatter of the population predictions as well as the large interindividual variability of the parameters suggests that ventilatory depression in individuals may be poorly predicted by population estimates.

The ability of opioids to provide profound analgesia is limited by the ventilatory depression associated with opioid overdose. Therefore, quantitation of opioid-induced ventilatory depression is a pharmacokinetic–pharmacodynamic problem relevant to the practice of anesthesia. In general, pharmacokinetic–pharmacodynamic modeling serves two purposes: identification of parameters to characterize the potency of drugs (typically, the concentration that causes half maximal effect, C⁵⁰) and identification of models that can be used for optimization of dosing regimens. The former can be achieved with steady state designs. However, the latter requires non–steady state data and, even in absence of feedback mechanisms, models to compensate for the hysteresis between the time courses of plasma concentration and effect. Unfortunately for the scientist, and fortunately for humans as biologic systems, ventilation is actively controlled by several feedback loops, including hypoxic ventilatory drive, hypercarbic ventilatory drive, and “wakefulness drive” (direct neural influences on respiratory activity). 13These sources of ventilatory drive contribute individually to the homeostasis of the respective variables (pH, Pao², Paco²). This poses additional difficulties for non–steady state models of drug-induced ventilatory depression. To model a feedback system, it must first be reduced to a manageable number of independent variables. Second, simple and stable relations of drugs on the feedback circuit must be characterized.

To determine the C⁵⁰of a drug influencing the behavior of a feedback-controlled system is a formidable task. There are two basic paradigms: the steady state approach and the dynamic or non–steady state approach. Maximal control of the experimental environment (clamped partial pressure of oxygen [Po²] and Pco²) and thereby opening the feedback loop collapses the problem to a simple, direct relation between dependent variable (MV) and independent variable (drug concentration), leading to a reliable estimate of the C⁵⁰in question obtainable with standard pharmacokinetic–pharmacodynamic models. The number of assumptions and the danger of model misspecification are negligible. Hence, the ventilatory depressant potency can be determined with a high degree of certainty, and precision and C⁵⁰values determined by this methodology are the “gold standard.” However, by eliminating the dynamic and feedback elements of the system, the approach cannot describe the time course of the effect in non–steady state situations, which includes the typical clinical applications of the model. The C⁵⁰refers to the concentration that decreases MV by 50% while maintaining a constant Paco², a situation never encountered in the clinical setting. An example of this steady state approach is the isohypercapnic approach to model drug-induced ventilatory depression, and this approach has successfully characterized the ventilatory depressant potency of fentanyl and morphine in neonatal dogs, 3alfentanil in humans, 4,6,7and remifentanil in humans. 4,5 

The alternative approach uses “real-world” dynamic data with non–steady state Paco²and drug concentration values and extracts the potency parameter from a composite model of the underlying system (the simplified  carbon dioxide kinetics, the simplified  carbon dioxide feedback loop, the pharmacokinetics of the drug, and the drug pharmacodynamics). With such a model, the number of assumptions and the danger of model misspecification are considerable. However, by including the feedback loop, the approach has the potential to describe the time course of the effect in non–steady state situations and is therefore useful for dose finding and clinical simulation and control. The system can be validated, at least in part, by comparing the potency parameter derived from non–steady state approaches with the potency estimated using steady state data, because the steady state model is a subset of the dynamic model.

This approach applied to drug-induced ventilatory depression has been called the modified indirect response model. It has been used successfully to determine the ventilatory depressant potency of alfentanil from non–steady state data, during and after a zero-order infusion (indirect response C⁵⁰60.3 ng/ml 6; isohypercapnic C⁵⁰s 49.5 and 75.3 ng/ml 4,7).

A dosing regimen designed for maximal “disturbance” of the system should be used in non–steady state studies to carefully assess the response. For this study, we used stepwise increases in remifentanil plasma concentrations using a target-controlled infusion device, up to an individually determined concentration associated with severe ventilatory depression. We expected the model to require a first-order equilibration delay to compensate for the time course of remifentanil transfer from the plasma to the site of remifentanil drug effect. The inclusion of a remifentanil plasma effect-site equilibration delay did not improve the quality of the fit. Our guess is that remifentanil's equilibration with the brain was so much faster than the equilibration of carbon dioxide that this additional delay was invisible to our measurements. For drugs with much slower equilibration between the plasma and the site of drug effect (e.g. , fentanyl, morphine), it would be critically important to base the modeling on effect-site opioid concentrations.

Because it is impossible to simultaneously determine carbon dioxide dynamics and pharmacodynamics from data with constantly changing drug and carbon dioxide concentrations over time, the experimental design included a drug-naive non–steady state carbon dioxide response curve in all subjects. We preferred Peco²to Paco²values for determination of the carbon dioxide response due to the higher resolution of the Peco²and the very small bias encountered (1.98–3.54 mmHg, 95% confidence interval).

Classically, carbon dioxide dynamics have been determined from the rebreathing design proposed by Read. 8In Read's approach, data analysis was limited to the linear portion of the stimulated hypercapnic carbon dioxide–versus –ventilation relation above the metabolic hyperbola. This approach has been shown to yield different carbon dioxide dynamics from the steady state design. 14We believe that the difference between classic Read rebreathing models of carbon dioxide dynamics and steady state carbon dioxide dynamics is at least partially due to an unaccounted hysteresis in the carbon dioxide–versus –ventilatory response relation. For this reason, we have used a modified rebreathing design in this study but added a first-order constant to permit modeling of the hysteresis between Paco²and ventilatory response. We would like to stress that our results are not intended to resolve the question about whether hysteresis explains discrepancies between the classic Read methodology and more contemporary steady state approaches to carbon dioxide dynamics. However, the good predictions of ventilation with non–steady state carbon dioxide obtained with our empirical and parsimonious model are consistent with the hypothesis that hysteresis is, at least partially, the explanation for the differences.

The portion of the carbon dioxide–versus –ventilatory response relation in the hypercapnic subject above approximately 45 mmHg is well approximated by a linear model. As the concentrations approach the metabolic hyperbola, the relation becomes flatter, with a distinct “hockey stick” appearance. Respiratory physiologists use two to three linear segments with corresponding thresholds to describe the effect of carbon dioxide on MV at or below resting values. 15,16Our nonlinear model is an attempt to create the same basic relation using a mathematically parsimonious approach. As mentioned before, we have not attempted to resolve this question with this pharmacodynamic study, but the good performance of the model suggests that our nonlinear equation is a reasonable approximation of the underlying relation. Interested readers are welcome to test different models for the shape of the carbon dioxide–versus –ventilation relation against a simulator that can be downloaded. **Our model was able to describe the observations in the simulated system well, although the F value we obtained from this simulator was approximately 6.

The core of our dynamic model of carbon dioxide is an indirect response model, which has been described previously. 6The carbon dioxide kinetic model and its parameters can be checked for plausibility by predicting the rate of increase of Paco²during apnea. The model predicts that Paco²increases 3.4 mmHg/min during apnea, which is in good agreement with the rate of increase during carbon dioxide rebreathing (3.7–4.2 mmHg/min [95% confidence interval]) and the standard view in respiratory physiology that carbon dioxide increases during apnea at 3–6 mmHg/min. 16The pharmacodynamic model can be checked for plausibility by comparing the C⁵⁰value with isohypercapnic C⁵⁰values. The C⁵⁰of 0.92 ng/ml is in good agreement with published isohypercapnic values of 1.12 ng/ml 5and 1.17 ng/ml. 4Thus, even though our carbon dioxide dynamic methodology differs from steady state approaches in several important details, our results are consistent with those obtained using steady state approaches.

The practical implications of our modeling exercise have been summarized in a simulation (fig. 4). Equal concentrations of remifentanil lead to different degrees of ventilatory depression, depending on the administration schedule. A 70-μg intravenous bolus of remifentanil causes a decrease of ventilation to essentially apneic values. The increase in Paco²occurs too slowly to adequately counteract the opioid effect. In fact, Paco²still increases after the nadir of MV has passed, demonstrating the inertia of the system. For an infusion designed to achieve the identical peak concentration (top panel ), the maximum depression of ventilation is much less because of the beneficial effect of the concurrent increase in Paco². Our model supports the clinical guidance that remifentanil boluses are inappropriate for patients when spontaneous respiration is desired (e.g. , for conscious sedation).

Finally, we would like to explain why patients experience only minimal impairment of steady state MV at the C⁵⁰for ventilatory depression and demonstrate that, under the assumption of constant carbon dioxide production, the steady state ventilation with unclamped carbon dioxide can directly be determined from the isohypercapnic/indirect response C⁵⁰value and the carbon dioxide sensitivity.

At steady state, with maintained spontaneous ventilation, the relation between Pecco²and (alveolar) ventilation is expressed by the metabolic hyperbola, regardless of the presence of drug. (Alveolar) ventilation equals the ratio of a constant (a) incorporating (constant) carbon dioxide production and Pecco²at the respective steady state and vice versa :

Combined with the equation for MV accounting for both carbon dioxide and drug effects (substituting Pecco²in equation 14into equation 11and solving for V̇^alv(ss)), we obtain:

which can be simplified to:

which permits us to express steady state ventilation as a fraction of baseline:

This equation can be used to determine the fractional steady state MV as a function of the drug concentration, isohypercapnic C⁵⁰, and F, the gain of the carbon dioxide–versus –ventilation response curve.

This equation also leads to an interesting observation about C⁵⁰. By definition, if the concentration of opioid equals the C⁵⁰, then ventilation will decrease by 50%, assuming no change in carbon dioxide. That is, of course, the acute response, not the steady state response. Based on the steady state equation above, when Cp(ss) = Cp⁵⁰, the equation reduces to:

Because F = 4.3, this can be solved for the fractional ventilation, which equals 88% of baseline. This explains why our dosing regimen, which was designed to avoid apnea, included the administration of concentrations fivefold higher than the C⁵⁰value. In fact, solving the equation for Cp(ss) = 5 × Cp⁵⁰yields a fractional alveolar ventilation of 67%, leading to a concomitant rise of fractional Paco²to 1.5 or, in other words, to 60 mmHg from a baseline of 40 mmHg. Of course, this control mechanism must collapse at higher drug concentrations or Paco²values, and we caution the reader not to extrapolate to higher than investigated concentration ranges. In addition, any decrease of the carbon dioxide production will lead to a further proportional decrease of fractional alveolar ventilation by shifting the position of the metabolic hyperbola to lower ventilation values.

We chose the commonly used sigmoid Emax model to relate opioid concentration to drug effect (ventilatory depression). Dahan et al.  7have typically modeled this relation using a power function. The power function performed relatively poorly describing the observations when compared to the sigmoid Emax model, demonstrating the importance of evaluating multiple models in pharmacodynamic research. However, neither model adequately captures periodic breathing and apnea in a meaningful way. The Emax model does not predict 0 ventilation at finite drug concentrations. The power model predicts ventilation up to fairly high concentrations, even though it is likely that patients intermittently become apneic at much lower concentrations. In our view, ventilation is a stochastic phenomenon, and apnea would be better predicted using a logistic probability model than a deterministic sigmoid Emax or power model.

In conclusion, we extended our previously described indirect response model to remifentanil in a dosing regimen consisting of multiple concentration steps. This indirect response model with carbon dioxide hysteresis yields estimates of C⁵⁰for ventilatory depression almost identical to those obtained with steady state methods. The dynamic carbon dioxide model may be useful in developing drug dosing regimens that minimize ventilatory depression, including the rational design of dosing regimens that balance the onset of opioid drug effect with ventilatory depression and rising carbon dioxide.