Background

Few studies address the dynamic effect of opioids on respiration. Models with intact feedback control of carbon dioxide on ventilation (non-steady-state models) that correctly incorporate the complex interaction among drug concentration, end-tidal partial pressure of carbon dioxide concentration, and ventilation yield reliable descriptions and predictions of the behavior of opioids. The authors measured the effect of remifentanil on respiration and developed a model of remifentanil-induced respiratory depression.

Methods

Ten male healthy volunteers received remifentanil infusions with different infusion speeds (target concentrations: 4-9 ng/ml; at infusion rates: 0.17-9 ng x ml x min) while awake and at the background of low-dose propofol. The data were analyzed with a nonlinear model consisting of two additive linear parts, one describing the depressant effect of remifentanil and the other describing the stimulatory effect of carbon dioxide on ventilation.

Results

The model adequately described the data including the occurrence of apnea. Most important model parameters were as follows: C50 for respiratory depression 1.6 +/- 0.03 ng/ml, gain of the respiratory controller (G) 0.42 - 0.1 l x min x Torr, and remifentanil blood effect site equilibration half-life (t(1/2)ke0) 0.53 +/- 0.2 min. Propofol caused a 20-50% reduction of C50 and G but had no effect on t(1/2)ke0. Apnea occurred during propofol infusion only. A simulation study revealed an increase in apnea duration at infusion speeds of 2.5-0.5 ng x ml x min followed by a reduction. At an infusion speed of < or = 0.31 ng x ml x min, no apnea was seen.

Conclusions

The effect of varying remifentanil infusions with and without a background of low-dose propofol on ventilation and end-tidal partial pressure of carbon dioxide concentration was described successfully using a non-steady-state model of the ventilatory control system. The model allows meaningful simulations and predictions.

  • ❖ The effects of opioids on respiration are opposed by stimulation of breathing from carbon dioxide retention, but most studies deal with detailed modeling of this relationship only at the steady state

  • ❖ In healthy human volunteers receiving remifentanil, respiration was modeled using a combination of depressant effects from the opioid and stimulation from carbon dioxide

  • ❖ Remifentanil infused in the presence of a fixed infusion of propofol resulted in apnea with increasing duration at remifentanil infusion speeds of 0.5–2.5 ng · ml−1· min−1

OPIOIDS affect breathing by activation of μ-opioid receptors expressed on respiratory neurons in the brainstem.1,2As a consequence, ventilation is depressed and arterial carbon dioxide increases. Because carbon dioxide activates chemoreceptors in the neck and brainstem (peripheral and central chemoreceptors), part of the opioid-induced respiratory depression is concealed by carbon dioxide-induced respiratory stimulation (the so-called carbon dioxide chemoreflex).3Especially, when the opioid slowly passes into the brain, the subsequent slow increase in carbon dioxide will offset major respiratory depression. Conversely, when the opioid rapidly passes the blood–brain barrier or the opioid is overdosed, depression of the respiratory neurons is faster and more noticeable than the respiratory stimulation from the carbon dioxide increase.3Then the effect of opioid is most dangerous. In general, opioids are considered safe with approximately 0.5% of patients receiving opioids for treatment of acute pain requiring immediate treatment for sometimes life-threatening respiratory depression.1However, in specific patient groups, this number is certainly much greater (e.g. , patients with sleep-related apnea, obesity, muscle weakness, and pulmonary disease). Furthermore, even in patients considered not at risk, opioid-induced mortality still occurs.4,5 

The number of studies on the effect of potent opioid analgesics on breathing is still limited. Even sparser are studies that quantify intravenous opioid effect on breathing using meaningfully parameterized pharmacodynamic models. The latter models are important as they allow, apart from the reliable description of opioid effect, the comparison among opioids (e.g. , on potency), the study of drug–drug interaction, and prediction of specific breathing-related idiosyncrasies such as the occurrence and duration of apnea. Current available models may be divided into two groups: (1) steady-state models, where the end-tidal partial pressure of carbon dioxide concentration (Pco2) is kept constant (by breath-to-breath manipulation of the inspired carbon dioxide concentration) and just the drug's effect on ventilation is measured (these models are also called open-loop models as the feedback loop between ventilation and arterial Pco2is broken—the loop is now open)6–10; and (2) non–steady-state models, in which the effects of the drug on arterial (or end-tidal) carbon dioxide and ventilation are both measured (these models are also called closed-loop models as the feedback loop between ventilation and arterial carbon dioxide remains active).11–13In contrast to steady-state models, non–steady-state models need to take into account the depressant effect that the opioid has on respiratory neurons in the brain causing the reduction of breathing and consequently the increase in arterial Pco2, but also and equally important, these models need to take into account the stimulatory effect of carbon dioxide on breathing. Only when both components are properly incorporated in the model, reliable estimates of the drug's respiratory potency are obtained, and a prediction of its respiratory behavior can be made.3 

We previously performed steady-state experiments and applied steady-state models to describe opioid-induced respiratory effects and their interaction with anesthetics (sevoflurane and propofol).9,10Although it allowed for the accurate description of the synergistic opioid-anesthetic interaction on breathing, this was unable to predict “real-life” non–steady-state conditions such as those occurring when drug concentrations rapidly change. In the current study, we performed non–steady-state experiments by applying increases in remifentanil concentration of different rates of rise. Experiments were performed in healthy volunteers in the awake condition and at the background of a low-dose propofol infusion. Next, we developed a non–steady-state pharmacokinetic–pharmacodynamic model of opioid-induced respiratory depression. The control of breathing is a complex system using both feedback and feedforward control tools to maintain cellular homeostasis.14Hence, it is important to make choices when considering the site of action of the opioid effect within the control system. We constructed a relatively simple model with drug concentration and end-tidal Pco2as input and measured inspired ventilation as output. Basic characteristics of the model are (1) it assumes that drug and carbon dioxide have opposing effects on breathing3; (2) it is based on the linear and well-described relationship between carbon dioxide and ventilation, V̇= G (Pco2− B), where V̇ is inspired minute ventilation, G is the gain of the ventilatory control system, and B the extrapolated end-tidal Pco2at which apnea occurs (apneic threshold)14–16; (3) the opioid effect is on B, whereas anesthetic effect is on G (fig. 1).10In our study, the model's behavior is tested to assess whether it accurately predicts apnea at finite opioid drug concentrations and whether all important model parameters are estimable (using a sensitivity analysis).

Fig. 1. (A ) Effect of remifentanil on the central gain of the respiratory controller as measured by Nieuwenhuijs et al.  10in steady-state experiments. At remifentanil concentrations more than 0, the gain remains constant at 1.6 ± 0.1 l · min−1· Torr−1. (B ) Effect of remifentanil on the apneic threshold (B) of the respiratory controller as measured by Nieuwenhuijs et al.  10B increases linearly with increasing remifentanil concentrations according to the function B(CREM) = 39 × (1 + CREM/6.14). (C ) Probability density function for end-tidal carbon dioxide pressure with ςA2 = 0.0225 and ςB2 = 0.25 (Eq. 8). CREM= remifentanil concentration in plasma.

Fig. 1. (A ) Effect of remifentanil on the central gain of the respiratory controller as measured by Nieuwenhuijs et al.  10in steady-state experiments. At remifentanil concentrations more than 0, the gain remains constant at 1.6 ± 0.1 l · min−1· Torr−1. (B ) Effect of remifentanil on the apneic threshold (B) of the respiratory controller as measured by Nieuwenhuijs et al.  10B increases linearly with increasing remifentanil concentrations according to the function B(CREM) = 39 × (1 + CREM/6.14). (C ) Probability density function for end-tidal carbon dioxide pressure with ςA2 = 0.0225 and ςB2 = 0.25 (Eq. 8). CREM= remifentanil concentration in plasma.

Close modal

Subjects

Ten healthy male volunteers (age, 18–30 yr; body mass index < 28 kg/m2) were recruited to participate in the study after approval of the protocol by the local Human Ethics Committee (Commissie Medische Ethiek, Leiden University Medical Center, Leiden, The Netherlands). Written and oral informed consent was obtained before inclusion in the study. All subjects were instructed not to eat or drink for at least 6 h before the study.

Study Design

After arrival in the laboratory, an arterial catheter for blood sampling was placed in the left or right radial artery. In the contralateral arm, an intravenous catheter was inserted for drug infusion. Each subject participated in three remifentanil infusions separated by 120-min washout intervals. The first two were without a background infusion of propofol; the last with a propofol infusion aimed at a target bispectral index (BIS) value of 80 (average target plasma concentration = 1 μg/ml). We randomly assigned one of the two initial experimental runs (i.e. , without propofol infusion) to be performed without blood sampling. In the other two runs, 3–6 arterial blood samples were obtained for remifentanil measurements at arbitrary time points.

Remifentanil was administered using a target-controlled infusion system. For remifentanil, we used a custom-built infusion pump that was programmed with a pharmacokinetic data set (Remifusor, University of Glasgow, Glasgow, United Kingdom).17We applied different remifentanil infusion schemes among the 10 subjects as described in table 1. We aimed at obtaining different rates of increases of the remifentanil plasma concentration (ranging from slow to fast). This was obtained by applying step increases in remifentanil plasma concentration (in 7 of 10 subjects, we used steps of 1 ng/ml) with varying step durations (0.5, 1, 1.5, 2, 3, 4, or 6 min) and with varying number of steps (2, 4, 5, or 6). In two subjects, we performed a single 1-min step with step sizes of 6 and 9 ng/ml. In the appropriate runs, one to three blood samples were randomly obtained during remifentanil infusion (but always just before a change in target remifentanil concentration); after infusion, two to three blood samples were obtained, again at random times. Each volunteer was subjected to three identical target remifentanil infusion schemes. In case of irregular breathing with periods of apnea (no breathing for periods > 10 s) and/or significant oxygen desaturations (< 90%), the subject was initially stimulated to take a deep breath. If this had no effect, the subject was artificially ventilated by a bag via  the mask and pneumotachograph for 20–30 s. The investigators could terminate or adapt the infusion at any time during the experiment when they felt that this was required to alleviate apnea and/or hypoxemia.

Table 1.  Remifentanil Infusion Schemes

Table 1.  Remifentanil Infusion Schemes
Table 1.  Remifentanil Infusion Schemes

Measurements

A facemask was applied over the mouth and nose. Inspired and expired gas flows were measured with a pneumotachograph connected to a pressure transducer and electronically integrated to yield a volume signal (Hans Rudolph, Myandotta, MI). During the studies, the subjects inhaled 100% oxygen. The oxygen and carbon dioxide concentrations of the inspired and expired gases and the arterial hemoglobin–oxygen saturation (Spo2) were measured with a Datex Multicap gas monitor (Datex-Engstrom, Helsinki, Finland). The electroencephalogram was recorded using an A-2000 monitor with software version 3.3 (Aspect Medical Systems, Norwood, MA). The monitor computed the BIS over 2-s epochs. We averaged the BIS values over 1 min. End-tidal Pco2and inspired minute ventilation were stored on a disc for further analysis. We further report the measured Spo2and BIS during the remifentanil infusions.

The samples for the determination of blood remifentanil concentrations were collected in tubes containing sodium heparin and immediately transferred to tubes containing 50% citric acid (to inactivate esterases) before freezing at −20°C. The assay method is based on tandem mass spectrometry detection.10 

Data Analysis

A population pharmacokinetic–pharmacodynamic analysis was performed on the data. The analysis was performed in two steps. In step 1, a population remifentanil pharmacokinetic analysis was performed. Next, by using the individual Bayesian pharmacokinetic estimates, a population pharmacodynamic analysis was performed.

Remifentanil Pharmacokinetics.

The description of remifentanil pharmacokinetics was aimed at obtaining individualized drug input functions to the pharmacologic model. The actual infusion rates from the log file of the target-controlled infusion device were used. The distributions of the structural parameters were fixed to the values of the three-compartmental model reported by Minto et al.  17(i.e. , the typical values and interindividual variabilities) and individualized by adjusting for age and lean body mass. A population analysis was performed, which allowed for Bayesian individualization of the structural parameters (albeit within the constrained distributions).18A remifentanil effect site was postulated, where the concentration lags with respect to the central compartment concentration as quantified by the equilibration half-life parameter t½ke0.

Carbon Dioxide Pharmacokinetics.

The relationship between carbon dioxide content (C) and its partial pressure (P) was assumed to be linear, so that P =λ0× C, where λ0= 0.863 Torr/(ml CO2in 100 ml blood) (fig. 2).19The following mass balance equations were used for the lungs and body (approximating the body by one compartment):

Fig. 2. Schematic representation of the pharmacokinetic–pharmacodynamic model. The model has three distinct parts. Part I is the remifentanil pharmacokinetic part, consisting of the distribution of remifentanil through the body, including the effect site (i.e.,  the respiratory controller in the brainstem) (via  a rate constant, ke0). Part II is the respiratory controller in the brainstem. Effect of remifentanil on the ventilatory control system results in a reduction of ventilation (via  a delay, τ). Part III is the part that describes carbon dioxide kinetics. Carbon dioxide production determines together with ventilation the arterial carbon dioxide concentration. Because remifentanil causes the reduction of ventilation and carbon dioxide production is minimally affected, the ability of the system to clear carbon dioxide has diminished, and arterial carbon dioxide concentrations increase. This will have a stimulatory effect on the respiratory controller (part II). The main effect of adding propofol on top of propofol is shown as an effect on the gain factor G. CL = clearance; CREM= remifentanil concentration in plasma; C50= concentration remifentanil causing 50% respiratory depression; G = gain of the ventilatory control system or slope of the hypercapnic ventilatory response; ke0= blood effect site rate constant; PE= end-tidal carbon dioxide partial pressure; PE_0= baseline (predrug) end-tidal carbon partial pressure; Q̇= cardiac output; τ= time constant of the ventilatory control system; V̇= inspired minute ventilation; 0 = baseline (predrug) inspired minute ventilation; CO2 = carbon dioxide production; VAL= alveolar volume; VTS= tissue volume; V1–3= volumes of compartments 1–3 of the kinetic remifentanil model.

Fig. 2. Schematic representation of the pharmacokinetic–pharmacodynamic model. The model has three distinct parts. Part I is the remifentanil pharmacokinetic part, consisting of the distribution of remifentanil through the body, including the effect site (i.e.,  the respiratory controller in the brainstem) (via  a rate constant, ke0). Part II is the respiratory controller in the brainstem. Effect of remifentanil on the ventilatory control system results in a reduction of ventilation (via  a delay, τ). Part III is the part that describes carbon dioxide kinetics. Carbon dioxide production determines together with ventilation the arterial carbon dioxide concentration. Because remifentanil causes the reduction of ventilation and carbon dioxide production is minimally affected, the ability of the system to clear carbon dioxide has diminished, and arterial carbon dioxide concentrations increase. This will have a stimulatory effect on the respiratory controller (part II). The main effect of adding propofol on top of propofol is shown as an effect on the gain factor G. CL = clearance; CREM= remifentanil concentration in plasma; C50= concentration remifentanil causing 50% respiratory depression; G = gain of the ventilatory control system or slope of the hypercapnic ventilatory response; ke0= blood effect site rate constant; PE= end-tidal carbon dioxide partial pressure; PE_0= baseline (predrug) end-tidal carbon partial pressure; Q̇= cardiac output; τ= time constant of the ventilatory control system; V̇= inspired minute ventilation; 0 = baseline (predrug) inspired minute ventilation; CO2 = carbon dioxide production; VAL= alveolar volume; VTS= tissue volume; V1–3= volumes of compartments 1–3 of the kinetic remifentanil model.

Close modal

where VALis the alveolar volume, PAis the arterial carbon dioxide pressure (which is assumed to equal alveolar pressure),   is the inspired minute ventilation,   is the cardiac output, PVis the venous carbon dioxide pressure, VTSis the apparent tissue volume, and CO2  is the carbon dioxide production. Because ventilation enters the model of carbon dioxide kinetics directly (see Eq. 1), no correction was made for dead space ventilation. Furthermore, λ1= k × PBW0/100 ≈ 10 and λ2= 100 ×λ0, where k is the volume conversion factor from standard temperature and pressure, dry to body temperature, and air saturated with water and PBWis the barometric pressure minus the pressure of air saturated with water. In the data analysis, we fixed VALto 3 l.19 CO2  was estimated from the baseline end-tidal carbon dioxide pressure and V̇. These baseline values,   and VTSwere parameters to be estimated.

Pharmacodynamic Analysis: Modeling the Effect of Remifentanil on the Ventilatory Control System.

The effect of end-tidal Pco2(PE) on ventilation under hyperoxic conditions can be modeled as follows (fig. 2)14–16:

where G is the central gain, B is the apneic threshold, and τ is a time constant. Nieuwenhuijs et al.  10characterized the remifentanil–propofol interaction on the ventilatory control system using a response surface modeling approach. For the current study, we reanalyzed those earlier data to characterize the effect of remifentanil on B and G (figs. 1A and B). From these analyses, we estimated that at remifentanil concentration (CREM) more than 0, G remained constant, whereas B increased linearly (the concentration remifentanil that doubles B = 6.14 ± 0.77 ng/ml). Hence, we assumed that in the current study, remifentanil changed B but not G.

In the steady state, we have:

with

where B0is the apneic threshold when CREM= 0 and C100the concentration remifentanil that causes a doubling of B. Rewriting Equation 3and defining baseline or resting end-tidal carbon dioxide concentration (i.e. , end-tidal Pco2before the remifentanil infusion) as PE_0, we get:

Baseline ventilation (0 ) = G · (PE_0− B0) and concentration remifentanil causing 50% respiratory depression, C50=

We then rewrite Equation 5into

Note that C50causes 50% depression of ventilation when G · (P − PE_0) = 0, which may occur when PE= PE_0(e.g. , after an acute remifentanil infusion when carbon dioxide did not rise as yet) or when G = 0 (as may occur when combining remifentanil with high propofol concentrations). This equation allows for the possibility of apnea at finite drug concentrations, which is appealing from a clinical point of view. Finally, in Equation 2, τ was fixed to 2.5 min.15 

Modeling the End-tidal Carbon Dioxide Pressure.

End-tidal Pco2is (under “normal” circumstances) an accurate indicator of alveolar Pco2. However, close to apnea, the breathing pattern is such that end-tidal Pco2as measured by the gas monitor is likely to be inaccurate and possibly measured too low. So, if we write for the residual error:

where PEis the measured value and E  the predicted value. The variance of ε should be smaller (ςA2) when PE> E  and larger when (ςB2) when PE< E . Therefore, the probability density of ε was written as:

which is a continuously differentiable function and integrates to 1. Because the distribution of ε is asymmetric and the mean of ε≠ 0, E  displayed in the figures is the mode. An example of the asymmetric probability density function with ςA2= 0.0225 and ςB2= 0.25 is given in figure 1C.

Furthermore, measured end-tidal carbon dioxide pressures were determined to be missing values if they were lower than 37.5 Torr or when the corresponding measured ventilation was less than 1 l/min because in those cases, it can be expected that end-tidal Pco2measurements are inaccurate.

Modeling Minute Ventilation.

Minute ventilation (V̇) was assumed to be normally distributed with variance ςv2. However, during apnea, manual ventilation (by mask) was applied. In that case, the measured ventilation values were determined to be missing values but used for the uptake and distribution model of carbon dioxide.

Statistical Analysis.

The models as described earlier were implemented in NONMEM VII (ICON Development Solutions, Ellicott City, MD).18The differential equations were solved with NONMEM's routine ADVAN6; the probability density functions (Eq. 8) were used with NONMEM'S LIKELIHOOD option. NONMEM VII's Markov Chain Monte Carlo Bayesian analysis method was used for parameter estimation. This method yields probability distributions of the model parameters from which means, standard errors, and 95% CI can be obtained. Uninformative priors were used for the interindividual variability terms for the pharmacodynamic analysis (in the pharmacokinetic analysis these were fixed); no priors were required for the structural parameters because of the highly informative data. An interoccasion variability term was incorporated for both parameters, 0  and PE_0(as their product is related to carbon dioxide production). Interindividual variability terms with large SE (larger than the estimate) were removed from the model. The burn-in samples were tested for convergence (all parameters and objective function more than 20 iterations, each 50 iterations apart; P < 0.05); 1,000 iterations were used to obtain parameter distributions. The significance of factors representing a deviance of pharmacokinetic parameters from those obtained by Minto et al.  17and significance of factors representing a decrease in the pharmacodynamic model parameters G and C50were tested by checking whether their 95% CIs excluded 1.

Sensitivity Analysis

A sensitivity analysis of a proposed model will indicate whether the parameter values of the model can be estimated with finite precision from the measured data.20The parameters may not be estimable for various reasons: because of the model structure, dependence on other parameters, or the specific input function chosen. We performed a sensitivity analysis of simulated data in which CREMincreases to 5 ng/ml in five steps using three distinct input functions: (1) step size = 1 ng/l, duration of step = 1 min; (2) step size = 1 ng/ml, duration of step = 5 min; and (3) step size = 1 ng/ml, duration of step = 0.1 min. The analysis was performed by fixing one parameter (i.e. , not allowing it to be estimated) at a time to a series of values (50–150%) of the “best” value of the parameter. Next, the other parameters were estimated, and the −2 log likelihood (−2LL) values were determined. This so-called likelihood profile method will show whether any of the parameters are estimable. If not, the curve of −2LL versus  the fixed parameter values (the “cost” function) will be flat.

Simulation Study

To get an indication on the effect of slowing the rate of rise of remifentanil on the nadir in ventilation or duration of apnea, we performed 30 simulations on remifentanil effect in the absence and presence of propofol. We simulated linear increases in remifentanil concentration at the effect site with rates of increase of 5–0.17 ng · ml−1· min−1. The infusion was terminated when the effect site concentration had reached 5 ng/ml.

All subjects completed the study without unintended effects. An example of one experimental run is given in figure 3; it is the data from one subject (id003) on the effect of a staircase increase in remifentanil concentration during a constant propofol infusion (run 3). Figure 3Ashows the target increase in remifentanil plasma concentration to 4 ng/ml. Note that the infusion was aborted early (because of the occurrence of apnea, fig. 3C). Figure 3Bshows the measured BIS values and fig. 3C, the inspired minute ventilation per breath. BIS values were on average 80 indicating moderate sedation, in agreement with the observation that the subject was unresponsive to verbal command. Breathing reduced rapidly on exposure to remifentanil, and apnea occurred after 3 min. Apnea was followed by irregular and cyclic breathing, which continued for 15–20 min, well after the remifentanil infusion was stopped. Note in fig. 3Dthat, at irregular breathing with low tidal volumes or apnea, an accurate measurement of end-tidal carbon dioxide was not possible.

Fig. 3. Example of the effect of remifentanil staircase infusion against the background of a constant propofol infusion. (A ) Target plasma remifentanil concentration (TCP). (B ) Bispectral index. (C ) Measured inspired ventilation (Vi). Note that ventilation quickly reaches apneic values after the initiation of the remifentanil infusion. Next, breathing remains irregular with reduced breathing frequency and periods of cyclic breathing. To get an indication of the sequential breathing pattern, the breaths are connected by a gray line . (D ) Measured end-tidal carbon dioxide concentration (PETCO2).

Fig. 3. Example of the effect of remifentanil staircase infusion against the background of a constant propofol infusion. (A ) Target plasma remifentanil concentration (TCP). (B ) Bispectral index. (C ) Measured inspired ventilation (Vi). Note that ventilation quickly reaches apneic values after the initiation of the remifentanil infusion. Next, breathing remains irregular with reduced breathing frequency and periods of cyclic breathing. To get an indication of the sequential breathing pattern, the breaths are connected by a gray line . (D ) Measured end-tidal carbon dioxide concentration (PETCO2).

Close modal

The infusion schemes applied in the studies are given in table 1. The estimated plasma concentration rates of increase varied from 0.17 to 9.0 ng · ml−1· min−1. BIS values were 93.2 ± 4.6 (mean ± SD) in remifentanil runs and 82.2 ± 4.8 (P < 0.05) in runs where remifentanil was given on top of propofol. The lowest values for Spo2were 93.8 ± 7.2% in the remifentanil runs and 87.5 ± 8.4% in the remifentanil–propofol runs (P < 0.05). Saturation values less than 90% occurred on average 30 ± 41 s in the remifentanil runs and 110 ± 86 s in the remifentanil–propofol runs (P = 0.05). Apnea did not occur in remifentanil runs but in eight remifentanil–propofol runs (averaged duration, 4.4 min; range, 1–7 min). In runs of subjects 007 and 009 (remifentanil rates of increase 0.17 and 0.22 ng · ml−1· min−1), the duration of apnea was at its lowest range, 1 and 2 min, respectively.

Pharmacokinetic Analysis

To get an indication of the goodness of the pharmacokinetic data fits, we plotted the measured concentrations versus  the individual and population-predicted remifentanil concentrations (figs. 4A and B). The plots indicate that the pharmacokinetic model adequately described the data. Two examples of pharmacokinetic data fits given in figure 4are in agreement with this statement. The pharmacokinetic parameter estimates were not significantly different from those of Minto et al.  17except for parameter V2(volume of compartment 2), which was a factor of 0.522 ± 0.125 (95% CI, 0.323–0.808) of that of Minto et al.  The factor in V2causes population measured versus  population-predicted concentrations to lie more closely on the line of identity. Without the factor, concentrations at the high end are underestimated and vice versa  (fig. 4C).

Fig. 4. Goodness-of-fit plots for the pharmacokinetic model. The measured concentrations (y-axis) versus  the individual- (A ) and population-predicted (B ) values. (C ) The pharmacokinetic analysis performed according to the parameter estimates of Minto et al.  17(i.e.  without a factor for V2).

Fig. 4. Goodness-of-fit plots for the pharmacokinetic model. The measured concentrations (y-axis) versus  the individual- (A ) and population-predicted (B ) values. (C ) The pharmacokinetic analysis performed according to the parameter estimates of Minto et al.  17(i.e.  without a factor for V2).

Close modal

Pharmacodynamic Analysis

The pharmacodynamic model adequately described the data. Examples of data fits are given in figure 5. The data are from one subject (id006) and are fits of ventilation and end-tidal carbon dioxide pressures under awake (fig. 5A–C) and sedated (fig. 5D–F) conditions. The effect of artificial ventilation by mask is clearly visible on end-tidal carbon dioxide in figure 5Fas these periods of artificial ventilation were incorporated in the pharmacodynamic model. Because no spontaneous ventilation was observable during the period of artificial ventilation, the fit through the ventilation data still (correctly) predicts apnea (fig. 5E).

Fig. 5. Examples of data fits of the effect of remifentanil on ventilation in one subject in the awake state (A –C ) and asleep with propofol (D –F ). (A, D ) Measured remifentanil concentration (closed circles ), pharmacokinetic data fits (thin line ), and estimated effect site (thick line ) concentrations. (B , E ) Minute ventilation with each spontaneous breath is depicted by open circles . Artificial breaths are depicted by closed squares . The lines  through the data are the data fits. During sleep (E ), the subject is apneic and requires artificial breathing assistance. (C , F ) End-tidal carbon dioxide values (open circles ) and data fits. (F ) The effect of the artificial ventilation is clearly visible in the data fits.

Fig. 5. Examples of data fits of the effect of remifentanil on ventilation in one subject in the awake state (A –C ) and asleep with propofol (D –F ). (A, D ) Measured remifentanil concentration (closed circles ), pharmacokinetic data fits (thin line ), and estimated effect site (thick line ) concentrations. (B , E ) Minute ventilation with each spontaneous breath is depicted by open circles . Artificial breaths are depicted by closed squares . The lines  through the data are the data fits. During sleep (E ), the subject is apneic and requires artificial breathing assistance. (C , F ) End-tidal carbon dioxide values (open circles ) and data fits. (F ) The effect of the artificial ventilation is clearly visible in the data fits.

Close modal

In table 2, the population pharmacodynamic model parameters are given. The concentration of remifentanil causing 50% respiratory depression is 1.6 ± 0.03 ng/ml. A notable observation is the low value for G in the remifentanil runs (0.42 l · min−1· Torr−1). Low-dose propofol significantly decreased parameters G by more than 50% to 0.19 l · min−1· Torr−1, and parameter C50by approximately 20% to 1.3 ng/ml. Propofol had no effect on baseline ventilation, baseline end-tidal carbon dioxide pressure, VTS, Q, and remifentanil t½ke0. All these latter values were within the expected ranges.

Table 2.  Pharmacodynamic Parameter Estimates

Table 2.  Pharmacodynamic Parameter Estimates
Table 2.  Pharmacodynamic Parameter Estimates

Sensitivity Analysis

The results of the likelihood profile method of the pharmacodynamic model are shown in figure 6. The Δ-2 log likelihood (or the cost function) indicates that the estimated model parameters (G, C50, t½ke0, VTS, Q) when applying relatively slow remifentanil input functions (that is, step durations of 1 and 5 min) were estimated with acceptable accuracy (±10% of the actual value). Because of noise on the simulated data, deviations from optimal parameter values were sometimes observed (i.e. , lower values of −2LL at “optimal” parameter values that were different from those used in the simulation). Much faster infusion (step duration is 0.1 min) yielded a significant estimation bias. Because we applied mostly slow input functions (step duration 1 min or larger in 9 of 10 subjects), we may assume that all important model parameters were estimable without any bias. Furthermore, visual inspection of the sensitivity analysis indicates that combining fast and slow input functions (as performed in our study) yields a reliable estimation without bias for all estimated parameters (intersection of all three in lines is around the x-value of 1, the simulated parameter value). Interestingly, Q was estimable at acceptable accuracy (which is an argument for the two-compartment carbon dioxide model that we used). Of the fixed parameters, estimation of parameters τ and VALVwas poor (τ) or impossible (VALV). We relate this to the specific design of the study. A different experiment design with periods of artificial ventilation will result in estimability of VALV, whereas steps in carbon dioxide would have resulted in the accurate estimation of τ.

Fig. 6. Sensitivity analysis for the model parameters. (A –D , F ) Estimated parameters and (E , G ) fixed parameters. Δ-2LL is the difference in −2 log likelihood between the simulated data and the best fit with one of the parameters held constant. The axis of each plot shows the parameter value as a fraction of the value observed or fixed in the analysis of the experimental data set. The short dashed line  indicates the P = 0.01 level, and the long dashed line  indicates the P = 0.05 level. •-• input function: step increase in plasma concentration of 1 ng/ml, duration of step is 1 min (rate of rise = 1 ng/ml per min), number of steps is 5. ○-○ input function: step increase in remifentanil concentration is 1 ng/ml, duration of step is 5 min (rate of rise = 0.2 ng/ml per min), number of steps is 5. x-x  input function: step increase in plasma concentration of 1 ng/ml, duration of step is 0.1 min (rate of rise = 10 ng/ml per min), number of steps is 5. C50= concentration remifentanil causing 50% respiratory depression; G = gain of the ventilatory control system or slope of the hypercapnic ventilatory response; ke0= blood effect site rate constant; τ= time constant of the ventilatory control system; Q = cardiac output; VALV= alveolar volume; VTS= tissue volume.

Fig. 6. Sensitivity analysis for the model parameters. (A –D , F ) Estimated parameters and (E , G ) fixed parameters. Δ-2LL is the difference in −2 log likelihood between the simulated data and the best fit with one of the parameters held constant. The axis of each plot shows the parameter value as a fraction of the value observed or fixed in the analysis of the experimental data set. The short dashed line  indicates the P = 0.01 level, and the long dashed line  indicates the P = 0.05 level. •-• input function: step increase in plasma concentration of 1 ng/ml, duration of step is 1 min (rate of rise = 1 ng/ml per min), number of steps is 5. ○-○ input function: step increase in remifentanil concentration is 1 ng/ml, duration of step is 5 min (rate of rise = 0.2 ng/ml per min), number of steps is 5. x-x  input function: step increase in plasma concentration of 1 ng/ml, duration of step is 0.1 min (rate of rise = 10 ng/ml per min), number of steps is 5. C50= concentration remifentanil causing 50% respiratory depression; G = gain of the ventilatory control system or slope of the hypercapnic ventilatory response; ke0= blood effect site rate constant; τ= time constant of the ventilatory control system; Q = cardiac output; VALV= alveolar volume; VTS= tissue volume.

Close modal

Simulation Study

The results of six simulations are given in figure 7. Linear remifentanil infusions with a rate of increase of 5 ng · ml−1· min−1(fig. 7A–C), 0.5 ng · ml−1· min−1(fig. 7D–F), and 0.2 ng · ml−1· min−1(fig. 7G–I) are shown for the awake condition (thin line) and at the background of propofol (thick lines). The effects of the rate of increase on the nadir in ventilation (in the awake studies) and duration of apnea (in the propofol studies) are given in figure 8for all 30 simulations. For both endpoints, the effect of slowing the rate of rise is biphasic. The nadir in ventilation decreased from 5 to 1 ng · ml−1· min−1(2.3–.6 l/min; fig. 5), after which it increased (5 ng · ml−1· min−1= linear infusion of 1 min; 1 ng · ml−1· min−1= linear infusion of 5 min). The duration of apnea increased from 2.5 to 0.6 ng · ml−1· min−1(2 min and 9 min, respectively) after which it decreased rapidly. At infusion rates of 0.31 ng · ml−1· min−1(16 min) and slower, no apnea occurred. During the two most rapid, short-term infusions (5 and 2.5 ng · ml−1· min−1or 1- and 2-min infusions), the remifentanil exposure was insufficient to cause apnea.

Fig. 7. Simulation study on the effect of changes in the rate of rise of the effect site remifentanil concentration on breathing (A , D , and G ) and end-tidal carbon dioxide pressure (B , E , and H ) in awake state (no propofol present, thin lines ) and during sleep (due to a low-dose propofol background infusion, thick line ). The remifentanil rates of rise are linear and vary from 5 ng · ml−1· min−1given for 1 min (A –C ) to 0.5 ng · ml−1· min−1given for 10 min (D –F ) and 0.2 ng · ml−1· min−1given for 25 min (G –I ), so that peak effect site remifentanil concentration was 5 ng/ml in all simulations. (C , F , and I ) The counterclockwise end-tidal carbon dioxide–ventilation loops (continuous lines ) and metabolic hyperbola (dashed line ). The gray triangles  depict the linear remifentanil infusion schemes.

Fig. 7. Simulation study on the effect of changes in the rate of rise of the effect site remifentanil concentration on breathing (A , D , and G ) and end-tidal carbon dioxide pressure (B , E , and H ) in awake state (no propofol present, thin lines ) and during sleep (due to a low-dose propofol background infusion, thick line ). The remifentanil rates of rise are linear and vary from 5 ng · ml−1· min−1given for 1 min (A –C ) to 0.5 ng · ml−1· min−1given for 10 min (D –F ) and 0.2 ng · ml−1· min−1given for 25 min (G –I ), so that peak effect site remifentanil concentration was 5 ng/ml in all simulations. (C , F , and I ) The counterclockwise end-tidal carbon dioxide–ventilation loops (continuous lines ) and metabolic hyperbola (dashed line ). The gray triangles  depict the linear remifentanil infusion schemes.

Close modal

Fig. 8. Simulation study on the effect of varying linear rates of rise of remifentanil effect site concentration on the nadir in ventilation (simulations performed in the absence of propofol, closed circles ) and on the duration of apnea (simulations performed in the presence of propofol, open squares ).

Fig. 8. Simulation study on the effect of varying linear rates of rise of remifentanil effect site concentration on the nadir in ventilation (simulations performed in the absence of propofol, closed circles ) and on the duration of apnea (simulations performed in the presence of propofol, open squares ).

Close modal

The carbon dioxide concentration–ventilation (counter clockwise) loops shown in figures 7C, F, and Iare graphical representations of the link between the two parameters in areas below and above the metabolic hyperbola (dashed lines). In the simulation studies, we assumed no effect of adding low dose propofol on metabolic rate (which was experimentally verified as resting Pco2and minute ventilation were similar between the awake and propofol sedates states). The graphs indicate that loops in the horizontal plane (such as observed in panel I) are desirable when aiming at and maintaining spontaneous breathing. The graphs show further that independent of the infusion rate hyperventilation in the recovery phase (upswing of the loop) is greater when ventilatory depression is more pronounced and carbon dioxide accumulates in the body.

Finally, the simulations are in close agreement with the observed data; for example, compare figure 7Bwith figure 5B, figure 7Ewith figure 5Eand figure 3C.

By using a “simple” non–steady-state pharmacokinetic– pharmacodynamic model of opioid-induced respiratory depression (Eq. 6), we described the ventilatory behavior of varying remifentanil infusion schemes in awake and propofol-sedated volunteers. The model is based on the linear relationship between carbon dioxide and ventilation. The most important parts of the model parameters were identifiable and estimable (e.g. , the drug concentration causing 50% respiratory depression, the gain factor of the respiratory controller, the remifentanil effect site equilibration half-life), whereas others were fixed to values obtained from previous studies from our laboratory or obtained from the literature (the time constant for carbon dioxide of the ventilatory control system and alveolar volume). As we assumed that ventilation was dependent not only on the remifentanil concentration at its effect site (the brainstem) but also on the metabolic product carbon dioxide, our model may be described as an indirect response model. The first (and only) previous indirect response model of opioid-induced respiratory depression was developed by Bouillon et al.  11,12Although their model differs at important points from ours, it even so shares important characteristics. We will discuss the similarities and differences between the two models in the section Model Comparisons later.

The Model: How It Works and Parameter Estimates

A schematic description of the model is given in figure 2. The model has three parts. The remifentanil pharmacokinetic part (part I), consists of the distribution of remifentanil throughout the body. Part of the remifentanil passes (with a delay described by rate constant ke0) to the effect site, the brainstem, where it affects the control of breathing (part II of the model via  term 1 of Equation 6: 0 [1 − CREM/2C50]). Respiration is diminished (with time constant τ) because of activation of μ-opioid receptors expressed on key parts of the ventilatory control system (e.g.,  premotor neurons of the ventral respiratory group [especially within the pre-Bötzinger complex] and the pontine respiratory group).1,2Part III of the model, the carbon dioxide kinetics, is affected by the diminished breathing as it reduces the efficacy of carbon dioxide output and as a consequence arterial Pco2increases. This again has an effect on the respiratory control system (part 2 of the model) as it stimulates breathing (via  term 2 of Equation 6: G ·[PE− PE_0]). So, two opposing additive effects influence breathing after remifentanil infusion: the direct depressant effect via  depression of respiratory neurons and a stimulatory effect of the increasing arterial Pco2. In figure 9A, the effect of just term 1 of Equation 6is plotted (lines 1 for awake and line 3 for asleep subjects). The effect of combining terms 1 and 2 is represented by lines 2 (awake) and 3 (asleep). It is apparent from the graphs that adding term 2 has a stimulatory effect on the remifentanil-ventilation data.

Fig. 9. (A ) The relationship between remifentanil concentration and ventilation is made up of two additive linear terms (Eq. 6). The term that describes the effect of remifentanil on ventilation when carbon dioxide has not accumulated (the first term of Eq. 6) is given by the linear line 1  for the awake state and line 3  for the propofol-sedated state. A steady state in ventilation occurs when carbon dioxide accumulates: for the awake state, this is reflected by the nonlinear line 2 , and for the sedated state by the nonlinear line 4 . Now both terms of Equation 6are active. Crosses  are C50values (C50is the [steady-state] concentration of remifentanil causing 50% respiratory depression). (B ) Comparison of the remifentanil–ventilation relationships derived from the model advanced by Bouillon et al.  12and our model. Lines 1  and 2  are the acute (line 1 ) and steady-state (line 2 ) relationships derived from Equation 6. Equivalent lines for the Bouillon model are lines 3  (acute relationship) and 4  (steady-state relationship) as derived from Equation 14 of Ref. 12. The large difference between models is the inability to predict apnea in the latter model (line 3 ). The steady-state relationships are alike (compare lines 2  and 4 ).

Fig. 9. (A ) The relationship between remifentanil concentration and ventilation is made up of two additive linear terms (Eq. 6). The term that describes the effect of remifentanil on ventilation when carbon dioxide has not accumulated (the first term of Eq. 6) is given by the linear line 1  for the awake state and line 3  for the propofol-sedated state. A steady state in ventilation occurs when carbon dioxide accumulates: for the awake state, this is reflected by the nonlinear line 2 , and for the sedated state by the nonlinear line 4 . Now both terms of Equation 6are active. Crosses  are C50values (C50is the [steady-state] concentration of remifentanil causing 50% respiratory depression). (B ) Comparison of the remifentanil–ventilation relationships derived from the model advanced by Bouillon et al.  12and our model. Lines 1  and 2  are the acute (line 1 ) and steady-state (line 2 ) relationships derived from Equation 6. Equivalent lines for the Bouillon model are lines 3  (acute relationship) and 4  (steady-state relationship) as derived from Equation 14 of Ref. 12. The large difference between models is the inability to predict apnea in the latter model (line 3 ). The steady-state relationships are alike (compare lines 2  and 4 ).

Close modal

An important assumption on which our model is based is that opioids (in our case remifentanil) cause a linear increase in the position of the ventilatory carbon dioxide response curve (in our model parameter B) with little change in the value of the response slope (in our model parameter G). There is ample evidence that opioids indeed cause a parallel shift of the steady-state V̇− Pco2response slope.10,21,22However, some studies indicate that there is some sex dependency with a reduction in slope in women (we exclusively performed studies in men), whereas others showed that the opioid effect is dependent on the technique used to measure the response slope.23As discussed previously, we consider the parallel shift of the V̇− Pco2response slope, a typical opioid effect and reduction of the slope because of a lowered arousal state (from sleep or sedatives/anesthetics).10 

The value of t½ke0that we observed (0.5 min) is in the same range as values from Babenco et al.  8using the isohypercapnic method (2 min) and studies on electroenecephalographic endpoints (1.6 min). The low t½ke0value is related to remifentanil's rapid passage across the blood–brain barrier. This rapid movement of remifentanil into the brain compartment is a potential danger because it may produce depression of respiratory neurons (via  term 1 of Eq. 6; fig. 9A) before any stimulatory effect of the accumulation of carbon dioxide is generated (via  term 2 of Eq. 6, fig. 9A). Other opioids with a much slower passage across the blood–brain barrier (such as morphine and its active metabolite morphine-6-glucuronide with a slow passage [t½ke0in the range of hours], but also drugs as buprenorphine [1 h], and fentanyl [5 min]) do allow for at least some accumulation of carbon dioxide while the respiratory neurons are being depressed, hence with a lesser chance of apnea.24–26However, this is only true when the drug is not overdosed. Otherwise, severe respiratory depression with apnea should always be in mind.4Our data further indicate that the chances of apnea are increased when the subject is asleep (BIS values = 80) with propofol.

When term 2 of Equation 6equals zero, the remifentanil concentration causing 50% depression of ventilation is by definition C50(now only term 1 of Equation 6is operative). This occurs in a situation in which Pco2remains constant at baseline levels (e.g. , in isohypercapnic experiments or when carbon dioxide has not yet risen above baseline values due to the rapid action of the opioid). Our C50(1.6 ng/ml) is, therefore, well comparable with values observed in isohypercapnic studies (compare with Babenco et al.  8who estimated a value of 1.4 ng/ml). In a previous isohypercapnic study from our laboratory studying the remifentanil propofol interaction, the C50value (0.7 ng/ml) for ventilation was measured at a fixed end-tidal Pco2of 55 Torr.10The observed C50is a low value, probably related to a differences in parameterization (see Eq. 4of Ref. 10) and the fact that we analyzed the whole remifentanil–propofol surface rather than just the remifentanil-ventilation relationship.

The value of G (table 2) is smaller than observed in a previous study on the effect of combining remifentanil and propofol on the ventilatory carbon dioxide response curve (value of G at CREM> 0 = 1.6 l · min−1· Torr−1; fig. 1).10The reason for the smaller estimated value of G may be the fact that we previously tested our subjects in normoxia versus  hyperoxia in the current study. Hyperoxia significantly blunts the peripheral chemoreceptors at the carotid bodies causing a 30–40% reduction of G.15Furthermore, we believe that the current study was performed on the dog-leg of the V̇− Pco2response slope. There is ample proof that the linear response curve flattens around resting carbon dioxide values.14The value observed by us is in agreement with the value presented by Babenco et al.  8after a remifentanil bolus infusion. Reassuring is the observation of a 50% reduction of G by low-dose propofol, which is in agreement with earlier findings.8 

The likelihood profile method is a tool to assess whether the parameters may be estimated from the data and also yield their 95% CIs. We observed that the ability to obtain accurate parameter estimates is dependent on the specific input function chosen. Slower remifentanil infusion rates (duration 1 min or longer) resulted in more precise estimates than a rapid bolus infusion (fig. 6). This may be related to the fact that we applied a rapid infusion in just two subjects, whereas the eight others received the slower rates. However, the analysis also indicates that when applying fast and slow infusion rates in one study and analyzing the data set using a population approach, the estimation precision is acceptable and that the fast infusions do not affect the outcome of the estimates negatively.

In summary, the estimated parameter values are in close agreement to previous reported values in the literature. This together with the results of the sensitivity analysis gives us confidence that our model adequately describes the effect of remifentanil on breathing. Furthermore, we were well able to describe and predict the occurrence of apnea.

Apnea

To the best of our knowledge, there are no earlier studies that assessed remifentanil effect on apnea in healthy volunteers. The study of Egan et al.  27comes closest to ours. They assessed the effect of various bolus doses of remifentanil and used a respiratory intervention scale, based primarily on Spo2, as endpoint. Low Spo2values (< 85%) did occur at the higher dose range (100 μg and greater) and predominantly in a population of older subjects (60–75 yr). The low Spo2values are most probably related to the occurrence of apnea. We simulated the dosing schedule of Egan et al.  and observed apnea at doses of 100 μg and greater (data not shown). The observation of hypoxemia is of interest. Similar to Egan et al. ,27we observed hypoxia although it was mild and occurred late in the experiment (3–5 min after the start of apnea). We relate this to the inspiration of 100% oxygen in our study, causing an oxygen store sufficient for 4–5 min before serious desaturation (Spo2< 90%) sets in. Low arterial oxygen has a stimulatory effect on breathing through activation of the peripheral chemoreceptors at the carotid bodies and increases the value of G.15,16This may have occurred in our experiments and may have had an effect on the duration of apnea. We did not take this into account in our current model but a G dependent on Spo2may be introduced in future models.

Remifentanil–Propofol Interaction

We previously assessed the interactive effects of remifentanil (0–2 ng/ml) and propofol (0–2 μg/ml) on breathing using steady-state isohypercapnic response surface modeling techniques.10We observed a synergistic interaction of the two drugs on resting ventilation, resting end-tidal Pco2, ventilation at a fixed Pco2of 55 Torr, and the V̇− Pco2response slope. Comparison with the current data analysis is difficult because of the differences in experimental set up and modeling approach. In contrast to our previous study, we currently studied just one low propofol plasma concentration (target on average 1 μg/ml). However, we believe (in retrospect) that the observation of large periods of apnea (1–7 min) during infusion of remifentanil against the background of low-dose propofol precludes testing of higher propofol doses in volunteers. The large difference in responses to remifentanil in awake and propofol-sedated states on apnea occurrence is in agreement with a synergistic interaction between the two drugs. Whether the propofol effect is related to γ-amino-butyric acidergic depression of respiratory neurons or secondary to the change in arousal-state remains unknown.2In agreement with the latter possibility is the finding that sleep induces a substantial enhancement of the depressant effect of morphine on ventilatory control.28 

Speed of Remifentanil Infusion

It has been suggested that by slowing the opioid infusion rate, the degree of respiratory depression diminishes because of the stimulatory effect of the accumulation of carbon dioxide.3,12We performed a simulation study to test this suggestion and observed a more complex interaction between remifentanil infusion rate and degree of depression of ventilation as measured by the nadir in ventilation (awake studies) and duration of apnea (propofol sedated studies; figs. 7 and 8). Going from a rapid (linear) infusion (5 ng · ml−1· min−1) to a slow (linear) infusion (0.17 ng · ml−1· min−1) in target effect site concentration (peak concentration = 5 ng/ml), both endpoints showed a worsening of ventilatory depression with a lower nadir in ventilation from 5 to 1 ng · ml−1· min−1(infusion duration: 1–5 min) and an increase in apnea duration from 2.5 to 0.55 ng · ml−1· min−1(infusion duration: 2–9 min). Thereafter, respiratory depression decreased with the slowing of the infusion rate. Apnea disappeared at infusion rates less than or equal to 0.31 ng · ml−1· min−1(infusion duration = 16 min). Our analysis indicates that the degree of respiratory depression (as defined by the nadir in ventilation and duration of apnea) is related to the opioid's pharmacokinetics, the speed of opioid infusion, the total amount of opioid given (at infusion durations of 1 and 2 min, the total amount of remifentanil given is insufficient to cause apnea; fig. 8), and the target plasma concentration, all relative to the carbon dioxide kinetics and dynamics. To prevent overt respiratory depression and apnea, one should consider all these variables. As the simulations performed by us are not to be extrapolated to other scenarios than those applied here, a simulation for each new circumstance should be performed.

Respiratory Variability

We observed great variability in the respiratory data during exposure to low-dose remifentanil (figs. 3 and 5). Mitsis et al.  29modeled variability of spontaneous respiration during low-dose remifentanil administration and concluded that the increase in variability because of the opioid was related to a decrease in the strength of the controller part of the ventilatory loop. Our observation of a low value of G is in agreement with this observation.

Model Comparisons

Bouillon et al.  11,12were the first to study and model the respiratory effects of opioids (alfentanil and remifentanil) in the non–steady state. Their initial attempts are well appreciated, and our modeling work should be considered as an extension to the original ideas postulated by Bouillon et al. 

Bouillon et al.  11–13used an indirect response model consisting of two multiplicative terms, a sigmoidal Emaxfunction and a power function of the form (Eq. 11 in ref. 12):

where F is the slope of the ventilatory carbon dioxide response curve (equivalent to our parameter G) and γ a shape parameter. We compared the model with our model in figure 9Bby plotting the acute (term 1, the sigmoid-Emaxfunction: line 3 in fig. 9B) and the combined (term 1 + term 2, line 4) steady-state relationships. The largest difference between the two models is the difference in the acute relationship as reflected by lines 1 and 3. Because of the sigmoidal Emaxnature of term 1 and the multiplicative nature of the model, no apnea may be described or predicted. At two times the C50, the Bouillon et al.  12model predicts ventilation levels more than 2 l/min, whereas in our approach, 2C50equals the concentration at which apnea occurs (CAPNEA= 3.2 ng/ml; fig. 9). The overall steady-state relationships (lines 2 and 4) are comparable between models. To get an impression of the ability of the model of Bouillon et al.  12to describe our data, we analyzed our data with their model. In contrast to Bouillon et al. , we added a delay between remifentanil concentration and the effect site (parameter ke0). Without this parameter, no meaningful data fits were obtained. Two main observations were made: as expected, the occurrence of apnea could not be modeled but yielded systematic misfits; and the effect of manual ventilation on Pco2yielded misfits as well. The latter is probably related to the single carbon dioxide compartment in the Bouillon et al.  model versus  two compartments in our model. In the study of Bouillon et al. , no systematic misfits were reported. However, this may partly be related to the remifentanil function applied (one to four 15-min steps of varying magnitudes of CREM, which yielded changes toward the steady state in both Pco2and ventilation) together with a relatively small number of arterial carbon dioxide samples. This may have yielded little contribution of the first term of the model to the overall effect but a predominant contribution of the power function, which depends on the accumulation of carbon dioxide. Possibly at different input functions (such as bolus infusions) or during sedation, the misspecifications of the model would become apparent. Bouillon et al.  12suggest to address the issue of apnea by using a logistic probability model. Our current model indicates that this is not required as apnea is accurately predicted at realistic drug concentrations. Our model has two linear terms (Eq. 6). To assess whether nonlinear functions would improve the data fits, we analyzed the data with a model consisting of power functions. No systematic improvements were observed (data not shown).

In conclusion, we developed a novel pharmacokinetic–pharmacodynamic model that is able to describe the non–steady-state effects of remifentanil on breathing under a variety of circumstances ranging from fast to slow drug infusions, under awake and sleeping conditions. Furthermore and possibly most important, the model allowed description and prediction of an important idiosyncrasy of the ventilatory control system, the development of opioid-induced apnea.

1.
Dahan A, Aarts L, Smith TW: Incidence, reversal, and prevention of opioid-induced respiratory depression. Anesthesiology 2010; 112:226–38
2.
Pattinson KTS: Opioids and the control of respiration. Br J Anaesth 2008; 100:747–58
3.
Dahan A, Teppema LJ: Influence of anaesthesia and analgesia on the control of breathing. Br J Anaesth 2003; 91:40–9
4.
Lötsch J, Dudziak R, Freynhagen R, Marschner J, Geisslinger G: Fatal respiratory depression after multiple intravenous morphine injections. Clin Pharmacokinet 2006; 45:1051–60
5.
Jumbelic MI: Deaths with transdermal fentanyl patches. Am J Forensic Med Pathol 2010; 31:1–4
6.
Blouin RT, Conrad PF, Gross JB: Time course of ventilatory depression following induction doses of propofol and thiopenthal. Anesthesiology 1991; 75:940–4
7.
Dershwitz M, Hoke FJ, Rosow CE, Michalowski P, Connors PM, Muir KT, Dienstag JL: Pharmacokinetics and pharmacodynamics of remifentanil in volunteer subjects with severe liver disease. Anesthesiology 1996; 84:812–20
8.
Babenco HD, Conrad PF, Gross JB: The pharmacodynamic effect of remifentanil bolus on ventilatory control. Anesthesiology 2000; 92:393–8
9.
Dahan A, Nieuwenhuijs D, Olofsen D, Sarton E, Romberg R, Teppema L: Response surface modeling of alfentanil-sevoflurane interaction on cardiorespiratory control and bispectral index. Anesthesiology 2000; 94:982–91
10.
Nieuwenhuijs D, Olofsen E, Romberg R, Sarton E, Ward DS, Engbers F, Vuyk J, Mooren R, Teppema L, Dahan A: Response surface modeling of remifentanil-propofol interaction on cardiorespiratory control and bispectral index. Anesthesiology 2003; 98:313–22
11.
Bouillon T, Schmidt C, Gartska G, Heimbach D, Stafforst D, Schwilden H, Hoeft A: Pharmacokinetic-pharmacodynamic modeling of the respiratory depressant effect of alfentanil. Anesthesiology 1999; 91:144–55
12.
Bouillon T, Bruhn J, Radu-Radulescu L, Andresen C, Cohane C, Shafer SL: A model of the ventilatory depressant potency of remifentanil in the non-steady state. Anesthesiology 2003; 99:779–87
13.
Bouillon T, Bruhn J, Radu-Radulescu L, Andresen C, Cohane C, Shafer SL: Mixed-effects modeling of the intrinsic ventilatory depressant potency of propofol in the non–steady state. Anesthesiology 2004; 100:240–50
14.
Cunningham DJC, Robbins PA, Wolff CB: Integration of respiratory responses to changes in alveolar partial pressures of CO2and O2and in arterial pH, Handbook of Physiology, Section 3: The Respiratory System, Volume II: Control of Breathing, Part 2. Edited by Fishman AP, Cherniack JG, Widdicombe JG, Geiger SR. Bethesda, American Physiological Society, 1986, pp 475–528Fishman AP, Cherniack JG, Widdicombe JG, Geiger SR
Bethesda
,
American Physiological Society
15.
Dahan A, DeGoede J, Berkenbosch A, Olievier I: The influence of oxygen on the ventilatory response to carbon dioxide in man. J Physiol (Lond) 1990; 428:485–99
16.
Dahan A, Nieuwenhuijs D, Teppema L: Plasticity of central chemoreceptors: Effect of bilateral carotid body resection on central CO2sensitivity. PLoS Med 2007; 4:e239
17.
Minto CF, Schnider TW, Egan TD, Youngs E, Lemmens HJM, Gambus PL, Billard V, Hoke JF, Moore KHP, Hermann DJ, Muir KT, Mandema JW, Shafer SL: Influence of age and gender on the pharmacokinetics and pharmacodynamics remifentanil: I. Model development. Anesthesiology 1997; 86:10–23
18.
Beal BL, Sheiner LB, Boeckman AJ: NONMEM User's Guide. Ellicott City, MD, Icon development Solutions, 1989–2006
19.
Lumb AB: Nunn's Applied Respiratory Physiology, 6th edition. Amsterdam, Elsevier Ltd., 2005
Amsterdam
,
Elsevier Ltd
20.
Ward DS, Dahan A, Mann CB: Modeling the dynamic ventilatory response to hypoxia in humans. Ann Biomed Eng 1992; 20:181–94
21.
Jordan C: Assessment of the effect of drugs on respiration. Br J Anaesth 1982; 54:763–82
22.
Bourke DL, Warley A: The steady-state and rebreathing methods compared during morfine administration in humans. J Physiol (Lond) 1989; 419:507–17
23.
Dahan A, Sarton E, Teppema L, Olievier C: Sex-related differences in the influence of morphine on ventilatory control in humans. Anesthesiology 1998; 88:903–13
24.
Yassen A, Olofsen E, Romberg R, Sarton E, Teppema L, Danhof M, Dahan A: Mechanism based PK/PD modeling of the respiratory depressant effect of buprenorphine and fentanyl in healthy volunteers. Clin Pharmacol Ther 2007; 81:50–8
25.
Sarton E, Olofsen E, Romberg R, den Hartigh J, Kest B, Nieuwenhuijs D, Burm A, Teppema L, Dahan A: Sex differences in morphine analgesia: An experimental study in healthy volunteers. Anesthesiology 2000; 95:1245–54
26.
Romberg R, Olofsen E, Sarton E, Teppema L, Dahan A: Pharmacodynamic effect of morphine-6-glucuronide versus  morphine on hypoxic and hypercapnic breathing in healthy volunteers. Anesthesiology 2003; 99:788–98
27.
Egan TD, Kern SE, Muir KT, White J: Remifentanil by bolus injector: A safety, pharmacokinetic, pharmacodynamic, and age effect investigation in human volunteers. Br J Anaesth 2004; 92:335–43
28.
Forrest WH, Bellville JW: The effect of sleep plus morphine on the respiratory response to carbon dioxide. Anesthesiology 1964; 25:137–41
29.
Mitsis GD, Governo RJM, Rogers R, Pattinson KTS: The effect of remifentanil on respiratory variability, evaluated with dynamic modeling. J Appl Physiol 2009; 106:1038–49