Free
Education  |   October 2003
A Model of the Ventilatory Depressant Potency of Remifentanil in the Non–steady State
Author Affiliations & Notes
  • Thomas Bouillon, M.D.
    *
  • Joergen Bruhn, M.D.
  • Lucian Radu-Radulescu, M.D.
  • Corina Andresen, M.D.
  • Carol Cohane, B.S.N.
    §
  • Steven L. Shafer, M.D.
  • *Staff Anesthesiologist, Department of Anesthesia, Inselspital Berne. †Assistant Professor, Department of Anesthesia, Bonn University Hospital, Bonn, Germany. ‡Research Associate, Department of Anesthesia, Duke University, Durham, North Carolina. §Clinical Research Coordinator, Department of Anesthesia, Palo Alto Health Care System and Stanford University School of Medicine. ∥Staff Anesthesiologist, Palo Alto Veterans Administration Health Care System; Professor of Anesthesia, Stanford University School of Medicine; Adjunct Professor of Biopharmaceutical Science, University of California at San Francisco, San Francisco, California.
  • Received from the Department of Anesthesia, Inselspital Berne, Berne, Switzerland; and the Palo Alto Veterans Administration Health Care System and Stanford University, Palo Alto, California.
Article Information
Education
Education   |   October 2003
A Model of the Ventilatory Depressant Potency of Remifentanil in the Non–steady State
Anesthesiology 10 2003, Vol.99, 779-787. doi:
Anesthesiology 10 2003, Vol.99, 779-787. doi:
This article is accompanied by an Editorial View. Please see: Gross JB: When you breathe in, you inspire; when you don't breathe, you expire: New insights regarding opioid-induced ventilatory depression. Anesthesiology 2003; 99:767–70.
THE ventilatory depressant action of μ-agonistic opioids represents their greatest potential source of toxicity. Application of morphine in the postoperative period resulted in significantly more frequent and severe decreases in oxygen saturation than regional anesthesia. 1 At a major Dutch University Hospital (Leiden, The Netherlands), 30% of all patients who received a general anesthetic for major surgery developed oxygen saturation below 90% in the postanesthesia care unit, caused at least partially by a blunted ventilatory drive. 2 As a result, there has been ongoing interest in finding methods that adequately describe opioid-induced ventilatory depression. Unfortunately, only on very few occasions were exact measures of potency (e.g.  , C50s) determined. 3–7 Most efforts to model the effects of opioids on ventilation have involved clamping partial pressure of carbon dioxide (Pco2) in a isohypercapnic approach, effectively canceling out carbon dioxide kinetics and dynamics and modeling drug-induced change of minute ventilation (MV) at constant Pco2. Because ventilatory depression in the clinical setting does not include carbon dioxide inhalation and carbon dioxide displays its own kinetics and dynamics, these models developed from isohypercapnic approaches cannot describe the extent and duration of opioid-induced ventilatory depression at non–steady state. A different approach accounting for both drug and carbon dioxide kinetics and dynamics has been suggested and successfully applied to the ventilatory depressant effect of alfentanil in the non–steady state. 6 The C50for alfentanil calculated from non–steady state data with an indirect response model (a model that does not relate [effect site] drug concentration to the value of a dependent variable but to the rate of change of a dependent variable) agreed with C50determinations from steady state approaches. 4,7 However, the non–steady state approach remains poorly tested (only a single zero-order infusion was used) and therefore requires further validation. We now apply the approach to non–steady state data for remifentanil, an opioid with rapid kinetics and effect compartment equilibration, which already has been investigated with isohypercapnic methods. 4,5 
Materials and Methods
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.
Subjects
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.
Study Design
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% CO2. 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).
Determination of Baseline Paco2
After a 5-min resting and equilibration period, an arterial blood sample for determination of arterial carbon dioxide tension (Paco2) was drawn.
Determination of the Ventilatory Response to Carbon Dioxide
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, 8 the 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 Paco2. After stabilization of MV at baseline levels, this part of the study was terminated.
Determination of Remifentanil-induced Ventilatory Depression
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 1running on a commercially available laptop computer. The remifentanil pharmacokinetic parameters were the covariate adjusted set reported by Minto et al.  9 Remifentanil was administered in ascending steps targeting plasma concentrations until end-tidal carbon dioxide (Peco2) 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 Paco2were drawn. Figure 1displays a typical example of this administration schedule and the corresponding Paco2values.
Table 1. Remifentanil Concentrations for the Ventilatory Depression Study
Image not available
Table 1. Remifentanil Concentrations for the Ventilatory Depression Study
×
Fig. 1. Experimental design displayed for one subject. After having undergone carbon dioxide rebreathing, the volunteer was subjected to stepwise increasing remifentanil concentrations applied with STANPUMP (3, 5, 7.5 ng/ml) targeting the plasma (solid line  ), not the effect compartment (dashed line  ). Unfilled circles  = remifentanil plasma concentration measurement; filled circles  = arterial carbon dioxide tension (PAco2) measurement.
Fig. 1. Experimental design displayed for one subject. After having undergone carbon dioxide rebreathing, the volunteer was subjected to stepwise increasing remifentanil concentrations applied with STANPUMP (3, 5, 7.5 ng/ml) targeting the plasma (solid line 
	), not the effect compartment (dashed line 
	). Unfilled circles 
	= remifentanil plasma concentration measurement; filled circles 
	= arterial carbon dioxide tension (PAco2) measurement.
Fig. 1. Experimental design displayed for one subject. After having undergone carbon dioxide rebreathing, the volunteer was subjected to stepwise increasing remifentanil concentrations applied with STANPUMP (3, 5, 7.5 ng/ml) targeting the plasma (solid line  ), not the effect compartment (dashed line  ). Unfilled circles  = remifentanil plasma concentration measurement; filled circles  = arterial carbon dioxide tension (PAco2) measurement.
×
Sampling and Data Processing
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 Paco2were 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 Paco2measurements were at best available every minute, MV measurements every 5 s and the further analysis required data pairs of Pco2and MV matched by time, Peco2values were used for calculation of the carbon dioxide response. To check for the validity of this approach, the average difference between Peco2and Paco2and the corresponding 95% confidence interval was calculated. Paco2measurements were used for calculation of drug-induced ventilatory depression.
Model Building
The pharmacokinetic model has been described previously. 10 In 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), V1= 13.50 l, V2= 8.64 l ± 27.4%, Cl1= 2.57 l/min ± 25.3%, Cl2= 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:MATHwhere θ(i)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 Paco2:MATHDVobsrefers to the observed value of the dependent variable (MV, V̇alv); DVexprefers 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.
Modeling the Carbon Dioxide Response Curves
Minute ventilation, MV, was modeled as a function of end-tidal carbon dioxide, Peco2, the independent variable. Hysteresis in the MV versus  Peco2relation was modeled using an effect compartment for carbon dioxide:MATHwhere Pecco2is the partial pressure of carbon dioxide at the effect site (“biophase”; mmHg) and ke0,CO2is the first-order equilibration constant between arterial and effect-site Pco2.
Before rebreathing, the system is at steady state, and baseline carbon dioxide, Pecco2(0), equals Peco2(0).
Although the relation between MV and Pecco2above the metabolic hyperbola (the steady state relation between alveolar ventilation and Paco2for a given carbon dioxide production and inspiratory fraction of carbon dioxide) can be well described by a straight line, 12,13 the 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, Pecco2was used as an independent variable of a nonlinear expression 6 :MATHwhere MV(Pecco2) is the MV depending on Pecco2(l/min); MV(0) is the baseline MV (l/min); Pecco2(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 Pecco2(t) and Pecco2(0). For reasons of completeness, a linear carbon dioxide response was also calculated as:MATHwhere SL is the slope of the carbon dioxide response curve. For each individual, baseline Pecco2was fixed to the measured value of baseline Peco2.
Modeling of Remifentanil-induced Ventilatory Depression
Because changes in Paco2during drug-induced ventilatory depression are not as fast as those during carbon dioxide rebreathing and Paco2is less artifact sensitive than Peco2(mask fit, shallow breathing), Paco2was used as the dependent variable for modeling of remifentanil-induced ventilatory depression. Although described previously, 6 the 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:MATHwhere Vdco2is the apparent volume of distribution of carbon dioxide (l), Paco2is the arterial partial pressure of carbon dioxide (mmHg), kinis the production rate of carbon dioxide (l/min), koutis the elimination rate of carbon dioxide (l/min), and 760 is the atmospheric pressure at sea level (mmHg). kout(t) can also be expressed as the product of alveolar ventilation (l/min) and the current Paco2divided by the barometric pressure, yielding:MATHUnder 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 Paco2and alveolar ventilation and becomes a constant:MATH
Rearranging for the change of Paco2, the dependent variable, over time yields:MATH
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:MATHwhere 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 ke0is 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 Pecco2on 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, cprefers to the plasma concentration of remifentanil, and C50refers to the remifentanil concentration at which V̇alvand therefore koutwill be decreased to 50% of the value in the absence of remifentanil, for unchanged Pecco2. 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 :
Results
General
One volunteer was excluded from the analysis because of a very low baseline Paco2, 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.
Carbon Dioxide Dynamics
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.
Table 2. Parameters of the Carbon Dioxide Response Curves Obtained in Absence of Drug, Nonlinear Carbon Dioxide Response
Image not available
Table 2. Parameters of the Carbon Dioxide Response Curves Obtained in Absence of Drug, Nonlinear Carbon Dioxide Response
×
Table 3. Parameters of the Carbon Dioxide Response Curves Obtained in Absence of Drug, Linear Carbon Dioxide Response
Image not available
Table 3. Parameters of the Carbon Dioxide Response Curves Obtained in Absence of Drug, Linear Carbon Dioxide Response
×
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.  Pco2in the effect compartment [Pecco2]). 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 ke0,CO2were used for the calculation of remifentanil-induced ventilatory depression.
Fig. 2. Carbon dioxide response curves. (Top  ) Bayesian predictions of minute ventilation versus  partial pressure of carbon dioxide (Pco2) in the effect compartment (n = 8). (Bottom  ) Goodness of fit for the nonlinear carbon dioxide response model. Predictions of minute ventilation based on typical values (unfilled circles  ) and Bayesian parameter estimates (filled circles  ) are plotted against measured minute ventilation. Note the pronounced interindividual variability in the carbon dioxide response. The Bayesian parameter values (individual carbon dioxide response) were used in the indirect response model.
Fig. 2. Carbon dioxide response curves. (Top 
	) Bayesian predictions of minute ventilation versus 
	partial pressure of carbon dioxide (Pco2) in the effect compartment (n = 8). (Bottom 
	) Goodness of fit for the nonlinear carbon dioxide response model. Predictions of minute ventilation based on typical values (unfilled circles 
	) and Bayesian parameter estimates (filled circles 
	) are plotted against measured minute ventilation. Note the pronounced interindividual variability in the carbon dioxide response. The Bayesian parameter values (individual carbon dioxide response) were used in the indirect response model.
Fig. 2. Carbon dioxide response curves. (Top  ) Bayesian predictions of minute ventilation versus  partial pressure of carbon dioxide (Pco2) in the effect compartment (n = 8). (Bottom  ) Goodness of fit for the nonlinear carbon dioxide response model. Predictions of minute ventilation based on typical values (unfilled circles  ) and Bayesian parameter estimates (filled circles  ) are plotted against measured minute ventilation. Note the pronounced interindividual variability in the carbon dioxide response. The Bayesian parameter values (individual carbon dioxide response) were used in the indirect response model.
×
Indirect Response Model Describing Drug-induced Ventilatory Depression
The remaining parameters characterizing the indirect response model describing remifentanil-induced ventilatory depression are summarized in table 4. Note that kel,CO2is the elimination constant of carbon dioxide and completely independent of drug action. Paco2(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.
Table 4. Pharmacodynamic Parameters of the Indirect Response Model
Image not available
Table 4. Pharmacodynamic Parameters of the Indirect Response Model
×
Figure 3displays the population and Bayesian prediction of Paco2for 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.
Fig. 3. (Top  ) Example of the fitted arterial carbon dioxide tension (Paco2) and concentration time course (same subject as shown in fig. 1). Dashed line  = plasma concentration as predicted by the target controlled infusion; unfilled circles  = prediction of plasma concentration based on Bayesian pharmacokinetic parameters (previously published pharmacokinetic model); filled circles  = measured Paco2; dotted line  = population prediction of Paco2; solid line  = Bayesian prediction of Paco2. (Bottom  ) Goodness of fit (all subjects). Predictions of Paco2based on typical values (unfilled circles  ) and Bayesian parameter estimates (filled circles  ) are plotted against measured Paco2.
Fig. 3. (Top 
	) Example of the fitted arterial carbon dioxide tension (Paco2) and concentration time course (same subject as shown in fig. 1). Dashed line 
	= plasma concentration as predicted by the target controlled infusion; unfilled circles 
	= prediction of plasma concentration based on Bayesian pharmacokinetic parameters (previously published pharmacokinetic model); filled circles 
	= measured Paco2; dotted line 
	= population prediction of Paco2; solid line 
	= Bayesian prediction of Paco2. (Bottom 
	) Goodness of fit (all subjects). Predictions of Paco2based on typical values (unfilled circles 
	) and Bayesian parameter estimates (filled circles 
	) are plotted against measured Paco2.
Fig. 3. (Top  ) Example of the fitted arterial carbon dioxide tension (Paco2) and concentration time course (same subject as shown in fig. 1). Dashed line  = plasma concentration as predicted by the target controlled infusion; unfilled circles  = prediction of plasma concentration based on Bayesian pharmacokinetic parameters (previously published pharmacokinetic model); filled circles  = measured Paco2; dotted line  = population prediction of Paco2; solid line  = Bayesian prediction of Paco2. (Bottom  ) Goodness of fit (all subjects). Predictions of Paco2based on typical values (unfilled circles  ) and Bayesian parameter estimates (filled circles  ) are plotted against measured Paco2.
×
Discussion
General
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, C50) 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). 13 These sources of ventilatory drive contribute individually to the homeostasis of the respective variables (pH, Pao2, Paco2). 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 C50of 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 [Po2] and Pco2) 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 C50in 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 C50values 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 C50refers to the concentration that decreases MV by 50% while maintaining a constant Paco2, 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, 3 alfentanil in humans, 4,6,7 and remifentanil in humans. 4,5 
The alternative approach uses “real-world” dynamic data with non–steady state Paco2and 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 C5060.3 ng/ml 6; isohypercapnic C50s 49.5 and 75.3 ng/ml 4,7).
Experimental and Model Building Considerations
Pharmacokinetics.
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.
Carbon Dioxide Dynamics.
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 Peco2to Paco2values for determination of the carbon dioxide response due to the higher resolution of the Peco2and 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. 8 In 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. 14 We 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 Paco2and 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,16 Our 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. 2Our model was able to describe the observations in the simulated system well, although the F value we obtained from this simulator was approximately 6.
Carbon Dioxide Kinetics and Pharmacodynamics.
The core of our dynamic model of carbon dioxide is an indirect response model, which has been described previously. 6 The carbon dioxide kinetic model and its parameters can be checked for plausibility by predicting the rate of increase of Paco2during apnea. The model predicts that Paco2increases 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. 16 The pharmacodynamic model can be checked for plausibility by comparing the C50value with isohypercapnic C50values. The C50of 0.92 ng/ml is in good agreement with published isohypercapnic values of 1.12 ng/ml 5 and 1.17 ng/ml. 4 Thus, 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 Paco2occurs too slowly to adequately counteract the opioid effect. In fact, Paco2still 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 Paco2. Our model supports the clinical guidance that remifentanil boluses are inappropriate for patients when spontaneous respiration is desired (e.g.  , for conscious sedation).
Fig. 4. The concentration–effect relation of remifentanil for ventilatory depression depends on the rate of application (simulation results). (Top  ) Concentration time course of remifentanil after a 70-μg bolus (solid line  ) and an infusion yielding a steady state concentration equal to the peak concentration after the bolus (dashed line  ). (Bottom  ) Corresponding time courses of arterial carbon dioxide tension (Paco2) and fractional alveolar ventilation. Solid line  = values corresponding to bolus application; dashed line  = values corresponding to infusion. The direction of change defines Paco2and ventilation plots (ventilation decreases, Paco2increases). Note that, because of the effect of carbon dioxide in the non–steady state, equal concentrations do not correspond to equal effect.
Fig. 4. The concentration–effect relation of remifentanil for ventilatory depression depends on the rate of application (simulation results). (Top 
	) Concentration time course of remifentanil after a 70-μg bolus (solid line 
	) and an infusion yielding a steady state concentration equal to the peak concentration after the bolus (dashed line 
	). (Bottom 
	) Corresponding time courses of arterial carbon dioxide tension (Paco2) and fractional alveolar ventilation. Solid line 
	= values corresponding to bolus application; dashed line 
	= values corresponding to infusion. The direction of change defines Paco2and ventilation plots (ventilation decreases, Paco2increases). Note that, because of the effect of carbon dioxide in the non–steady state, equal concentrations do not correspond to equal effect.
Fig. 4. The concentration–effect relation of remifentanil for ventilatory depression depends on the rate of application (simulation results). (Top  ) Concentration time course of remifentanil after a 70-μg bolus (solid line  ) and an infusion yielding a steady state concentration equal to the peak concentration after the bolus (dashed line  ). (Bottom  ) Corresponding time courses of arterial carbon dioxide tension (Paco2) and fractional alveolar ventilation. Solid line  = values corresponding to bolus application; dashed line  = values corresponding to infusion. The direction of change defines Paco2and ventilation plots (ventilation decreases, Paco2increases). Note that, because of the effect of carbon dioxide in the non–steady state, equal concentrations do not correspond to equal effect.
×
Finally, we would like to explain why patients experience only minimal impairment of steady state MV at the C50for 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 C50value and the carbon dioxide sensitivity.
At steady state, with maintained spontaneous ventilation, the relation between Pecco2and (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 Pecco2at the respective steady state and vice versa  :
Combined with the equation for MV accounting for both carbon dioxide and drug effects (substituting Pecco2in equation 14into equation 11and solving for V̇alv(ss)), we obtain:MATHwhich can be simplified to:MATHwhich permits us to express steady state ventilation as a fraction of baseline:MATH
This equation can be used to determine the fractional steady state MV as a function of the drug concentration, isohypercapnic C50, and F, the gain of the carbon dioxide–versus  –ventilation response curve.
This equation also leads to an interesting observation about C50. By definition, if the concentration of opioid equals the C50, 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) = Cp50, the equation reduces to:MATH
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 C50value. In fact, solving the equation for Cp(ss) = 5 × Cp50yields a fractional alveolar ventilation of 67%, leading to a concomitant rise of fractional Paco2to 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 Paco2values, 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.  7 have 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 C50for 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.
References
Catley DM, Thornton C, Jordan C, Lehane JR, Royston D, Jones JG: Pronounced, episodic oxygen desaturation in the postoperative period: Its association with ventilatory pattern and analgesic regimen. A nesthesiology 1985; 63: 20–8Catley, DM Thornton, C Jordan, C Lehane, JR Royston, D Jones, JG
Sarton E, Romberg R, Nieuwenhuijs D, Teppema LJ, Vuyk J, Schraag S, Dahan A: Einfluss von Anasthetika auf die Atemkontrolle [The effect of anesthetics on control of respiration]. Anaesthesist 2002; 51: 285–91Sarton, E Romberg, R Nieuwenhuijs, D Teppema, LJ Vuyk, J Schraag, S Dahan, A
Bragg P, Zwass MS, Lau M, Fisher DM: Opioid pharmacodynamics in neonatal dogs: Differences between morphine and fentanyl. J Appl Physiol 1995; 79: 1519–24Bragg, P Zwass, MS Lau, M Fisher, DM
Glass PS, Iselin-Chaves IA, Goodman D, Delong E, Hermann DJ: Determination of the potency of remifentanil compared with alfentanil using ventilatory depression as the measure of opioid effect. A nesthesiology 1999; 90: 1556–63Glass, PS Iselin-Chaves, IA Goodman, D Delong, E Hermann, DJ
Babenco HD, Conard PF, Gross JB: The pharmacodynamic effect of a remifentanil bolus on ventilatory control. A nesthesiology 2000; 92: 393–8Babenco, HD Conard, PF Gross, JB
Bouillon T, Schmidt C, Garstka G, Heimbach D, Stafforst D, Schwilden H, Hoeft A: Pharmacokinetic–pharmacodynamic modeling of the respiratory depressant effect of alfentanil. A nesthesiology 1999; 91: 144–55Bouillon, T Schmidt, C Garstka, G Heimbach, D Stafforst, D Schwilden, H Hoeft, A
Dahan A, Nieuwenhuijs D, Olofsen E, Sarton E, Romberg R, Teppema L: Response surface modeling of alfentanil–sevoflurane interaction on cardiorespiratory control and Bispectral Index. A nesthesiology 2001; 94: 982–91Dahan, A Nieuwenhuijs, D Olofsen, E Sarton, E Romberg, R Teppema, L
Read DJ: A clinical method for assessing the ventilatory response to carbon dioxide. Australas Ann Med 1967; 16: 20–32Read, DJ
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. A nesthesiology 1997; 86: 10–23Minto, 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
Bouillon T, Bruhn J, Radu-Radulescu L, Bertaccini E, Park S, Shafer SL: Non–steady state analysis of the pharmacokinetic interaction between propofol and remifentanil. A nesthesiology 2002; 97: 1350–62Bouillon, T Bruhn, J Radu-Radulescu, L Bertaccini, E Park, S Shafer, SL
Beal SL, Sheiner LB: NONMEM User's Guides. San Francisco, NONMEM Project Group, University of California, 1992
Dahan A, Berkenbosch A, DeGoede J, Olievier IC, Bovill JG: On a pseudo-rebreathing technique to assess the ventilatory sensitivity to carbon dioxide in man. J Physiol 1990; 423: 615–29Dahan, A Berkenbosch, A DeGoede, J Olievier, IC Bovill, JG
Longobardo G, Evangelisti CJ, Cherniack NS: Effects of neural drives on breathing in the awake state in humans. Respir Physiol 2002; 129: 317–33Longobardo, G Evangelisti, CJ Cherniack, NS
Berkenbosch A, Bovill JG, Dahan A, DeGoede J, Olievier IC: The ventilatory CO2sensitivities from Read's rebreathing method and the steady-state method are not equal in man. J Physiol 1989; 411: 367–77Berkenbosch, A Bovill, JG Dahan, A DeGoede, J Olievier, IC
Duffin J, Mohan RM, Vasiliou P, Stephenson R, Mahamed S: A model of the chemoreflex control of breathing in humans: Model parameters measurement. Respir Physiol 2000; 120: 13–26Duffin, J Mohan, RM Vasiliou, P Stephenson, R Mahamed, S
Nunn JF: Applied Respiratory Physiology, 4th edition. Oxford, England, Butterworth Heinemann, 1993, p 237
Fig. 1. Experimental design displayed for one subject. After having undergone carbon dioxide rebreathing, the volunteer was subjected to stepwise increasing remifentanil concentrations applied with STANPUMP (3, 5, 7.5 ng/ml) targeting the plasma (solid line  ), not the effect compartment (dashed line  ). Unfilled circles  = remifentanil plasma concentration measurement; filled circles  = arterial carbon dioxide tension (PAco2) measurement.
Fig. 1. Experimental design displayed for one subject. After having undergone carbon dioxide rebreathing, the volunteer was subjected to stepwise increasing remifentanil concentrations applied with STANPUMP (3, 5, 7.5 ng/ml) targeting the plasma (solid line 
	), not the effect compartment (dashed line 
	). Unfilled circles 
	= remifentanil plasma concentration measurement; filled circles 
	= arterial carbon dioxide tension (PAco2) measurement.
Fig. 1. Experimental design displayed for one subject. After having undergone carbon dioxide rebreathing, the volunteer was subjected to stepwise increasing remifentanil concentrations applied with STANPUMP (3, 5, 7.5 ng/ml) targeting the plasma (solid line  ), not the effect compartment (dashed line  ). Unfilled circles  = remifentanil plasma concentration measurement; filled circles  = arterial carbon dioxide tension (PAco2) measurement.
×
Fig. 2. Carbon dioxide response curves. (Top  ) Bayesian predictions of minute ventilation versus  partial pressure of carbon dioxide (Pco2) in the effect compartment (n = 8). (Bottom  ) Goodness of fit for the nonlinear carbon dioxide response model. Predictions of minute ventilation based on typical values (unfilled circles  ) and Bayesian parameter estimates (filled circles  ) are plotted against measured minute ventilation. Note the pronounced interindividual variability in the carbon dioxide response. The Bayesian parameter values (individual carbon dioxide response) were used in the indirect response model.
Fig. 2. Carbon dioxide response curves. (Top 
	) Bayesian predictions of minute ventilation versus 
	partial pressure of carbon dioxide (Pco2) in the effect compartment (n = 8). (Bottom 
	) Goodness of fit for the nonlinear carbon dioxide response model. Predictions of minute ventilation based on typical values (unfilled circles 
	) and Bayesian parameter estimates (filled circles 
	) are plotted against measured minute ventilation. Note the pronounced interindividual variability in the carbon dioxide response. The Bayesian parameter values (individual carbon dioxide response) were used in the indirect response model.
Fig. 2. Carbon dioxide response curves. (Top  ) Bayesian predictions of minute ventilation versus  partial pressure of carbon dioxide (Pco2) in the effect compartment (n = 8). (Bottom  ) Goodness of fit for the nonlinear carbon dioxide response model. Predictions of minute ventilation based on typical values (unfilled circles  ) and Bayesian parameter estimates (filled circles  ) are plotted against measured minute ventilation. Note the pronounced interindividual variability in the carbon dioxide response. The Bayesian parameter values (individual carbon dioxide response) were used in the indirect response model.
×
Fig. 3. (Top  ) Example of the fitted arterial carbon dioxide tension (Paco2) and concentration time course (same subject as shown in fig. 1). Dashed line  = plasma concentration as predicted by the target controlled infusion; unfilled circles  = prediction of plasma concentration based on Bayesian pharmacokinetic parameters (previously published pharmacokinetic model); filled circles  = measured Paco2; dotted line  = population prediction of Paco2; solid line  = Bayesian prediction of Paco2. (Bottom  ) Goodness of fit (all subjects). Predictions of Paco2based on typical values (unfilled circles  ) and Bayesian parameter estimates (filled circles  ) are plotted against measured Paco2.
Fig. 3. (Top 
	) Example of the fitted arterial carbon dioxide tension (Paco2) and concentration time course (same subject as shown in fig. 1). Dashed line 
	= plasma concentration as predicted by the target controlled infusion; unfilled circles 
	= prediction of plasma concentration based on Bayesian pharmacokinetic parameters (previously published pharmacokinetic model); filled circles 
	= measured Paco2; dotted line 
	= population prediction of Paco2; solid line 
	= Bayesian prediction of Paco2. (Bottom 
	) Goodness of fit (all subjects). Predictions of Paco2based on typical values (unfilled circles 
	) and Bayesian parameter estimates (filled circles 
	) are plotted against measured Paco2.
Fig. 3. (Top  ) Example of the fitted arterial carbon dioxide tension (Paco2) and concentration time course (same subject as shown in fig. 1). Dashed line  = plasma concentration as predicted by the target controlled infusion; unfilled circles  = prediction of plasma concentration based on Bayesian pharmacokinetic parameters (previously published pharmacokinetic model); filled circles  = measured Paco2; dotted line  = population prediction of Paco2; solid line  = Bayesian prediction of Paco2. (Bottom  ) Goodness of fit (all subjects). Predictions of Paco2based on typical values (unfilled circles  ) and Bayesian parameter estimates (filled circles  ) are plotted against measured Paco2.
×
Fig. 4. The concentration–effect relation of remifentanil for ventilatory depression depends on the rate of application (simulation results). (Top  ) Concentration time course of remifentanil after a 70-μg bolus (solid line  ) and an infusion yielding a steady state concentration equal to the peak concentration after the bolus (dashed line  ). (Bottom  ) Corresponding time courses of arterial carbon dioxide tension (Paco2) and fractional alveolar ventilation. Solid line  = values corresponding to bolus application; dashed line  = values corresponding to infusion. The direction of change defines Paco2and ventilation plots (ventilation decreases, Paco2increases). Note that, because of the effect of carbon dioxide in the non–steady state, equal concentrations do not correspond to equal effect.
Fig. 4. The concentration–effect relation of remifentanil for ventilatory depression depends on the rate of application (simulation results). (Top 
	) Concentration time course of remifentanil after a 70-μg bolus (solid line 
	) and an infusion yielding a steady state concentration equal to the peak concentration after the bolus (dashed line 
	). (Bottom 
	) Corresponding time courses of arterial carbon dioxide tension (Paco2) and fractional alveolar ventilation. Solid line 
	= values corresponding to bolus application; dashed line 
	= values corresponding to infusion. The direction of change defines Paco2and ventilation plots (ventilation decreases, Paco2increases). Note that, because of the effect of carbon dioxide in the non–steady state, equal concentrations do not correspond to equal effect.
Fig. 4. The concentration–effect relation of remifentanil for ventilatory depression depends on the rate of application (simulation results). (Top  ) Concentration time course of remifentanil after a 70-μg bolus (solid line  ) and an infusion yielding a steady state concentration equal to the peak concentration after the bolus (dashed line  ). (Bottom  ) Corresponding time courses of arterial carbon dioxide tension (Paco2) and fractional alveolar ventilation. Solid line  = values corresponding to bolus application; dashed line  = values corresponding to infusion. The direction of change defines Paco2and ventilation plots (ventilation decreases, Paco2increases). Note that, because of the effect of carbon dioxide in the non–steady state, equal concentrations do not correspond to equal effect.
×
Table 1. Remifentanil Concentrations for the Ventilatory Depression Study
Image not available
Table 1. Remifentanil Concentrations for the Ventilatory Depression Study
×
Table 2. Parameters of the Carbon Dioxide Response Curves Obtained in Absence of Drug, Nonlinear Carbon Dioxide Response
Image not available
Table 2. Parameters of the Carbon Dioxide Response Curves Obtained in Absence of Drug, Nonlinear Carbon Dioxide Response
×
Table 3. Parameters of the Carbon Dioxide Response Curves Obtained in Absence of Drug, Linear Carbon Dioxide Response
Image not available
Table 3. Parameters of the Carbon Dioxide Response Curves Obtained in Absence of Drug, Linear Carbon Dioxide Response
×
Table 4. Pharmacodynamic Parameters of the Indirect Response Model
Image not available
Table 4. Pharmacodynamic Parameters of the Indirect Response Model
×