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.

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.

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.

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.

## What We Already Know about This Topic

❖ 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

## What This Article Tells Us That Is New

❖ 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 (Pco^{2}) 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 Pco^{2}is 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 Pco^{2}, 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 Pco^{2}as 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 (Pco^{2}− B), where V̇ is inspired minute ventilation, G is the gain of the ventilatory control system, and B the extrapolated end-tidal Pco^{2}at 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).

## Materials and Methods

### Subjects

Ten healthy male volunteers (age, 18–30 yr; body mass index < 28 kg/m^{2}) 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.

### 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 (Spo^{2}) 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 Pco^{2}and inspired minute ventilation were stored on a disc for further analysis. We further report the measured Spo^{2}and 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^{½}k^{e0}.

#### 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 CO^{2}in 100 ml blood) (fig. 2).19The following mass balance equations were used for the lungs and body (approximating the body by one compartment):

where V^{AL}is the alveolar volume, P^{A}is the arterial carbon dioxide pressure (which is assumed to equal alveolar pressure), **V̇** is the inspired minute ventilation, **Q̇** is the cardiac output, P^{V}is the venous carbon dioxide pressure, V^{TS}is the apparent tissue volume, and **V̇ ^{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 × P

^{BW}/λ

^{0}/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 P

^{BW}is the barometric pressure minus the pressure of air saturated with water. In the data analysis, we fixed V

^{AL}to 3 l.19

**V̇**was estimated from the baseline end-tidal carbon dioxide pressure and V̇. These baseline values,

^{CO2}**Q̇**and V

^{TS}were parameters to be estimated.

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

The effect of end-tidal Pco^{2}(P^{E}) 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 (C^{REM}) 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 B^{0}is the apneic threshold when C^{REM}= 0 and C^{100}the concentration remifentanil that causes a doubling of B. Rewriting Equation 3and defining baseline or resting end-tidal carbon dioxide concentration (*i.e.* , end-tidal Pco^{2}before the remifentanil infusion) as P^{E_0}, we get:

Baseline ventilation (**V̇ ^{0}** ) = G · (P

^{E_0}− B

^{0}) and concentration remifentanil causing 50% respiratory depression, C

^{50}=

We then rewrite Equation 5into

Note that C^{50}causes 50% depression of ventilation when G · (P − P^{E_0}) = 0, which may occur when P^{E}= P^{E_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 Pco^{2}is (under “normal” circumstances) an accurate indicator of alveolar Pco^{2}. However, close to apnea, the breathing pattern is such that end-tidal Pco^{2}as measured by the gas monitor is likely to be inaccurate and possibly measured too low. So, if we write for the residual error:

where P^{E}is the measured value and **P̂ ^{E}** the predicted value. The variance of ε should be smaller (ς

^{A}

^{2}) when P

^{E}>

**P̂**and larger when (ς

^{E}^{B}

^{2}) when P

^{E}<

**P̂**. Therefore, the probability density of ε was written as:

^{E}which is a continuously differentiable function and integrates to 1. Because the distribution of ε is asymmetric and the mean of ε≠ 0, **P̂ ^{E}** displayed in the figures is the mode. An example of the asymmetric probability density function with ς

^{A}

^{2}= 0.0225 and ς

^{B}

^{2}= 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 Pco^{2}measurements are inaccurate.

#### Modeling Minute Ventilation.

Minute ventilation (V̇) was assumed to be normally distributed with variance ς^{v}^{2}. 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, **V̇ ^{0}** and P

^{E_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 C

^{50}were 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 C^{REM}increases 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.

## Results

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.

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 Spo^{2}were 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 V^{2}(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 V^{2}causes 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).

### 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).

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 C^{50}by approximately 20% to 1.3 ng/ml. Propofol had no effect on baseline ventilation, baseline end-tidal carbon dioxide pressure, V^{TS}, Q, and remifentanil t^{½}k^{e0}. All these latter values were within the expected ranges.

### 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, C^{50}, t½k^{e0}, V^{TS}, 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 V^{ALV}was poor (τ) or impossible (V^{ALV}). We relate this to the specific design of the study. A different experiment design with periods of artificial ventilation will result in estimability of V^{ALV}, whereas steps in carbon dioxide would have resulted in the accurate estimation of τ.

### 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^{−1}or 1- and 2-min infusions), the remifentanil exposure was insufficient to cause apnea.

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 Pco^{2}and 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.

## Discussion

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 k^{e0}) to the effect site, the brainstem, where it affects the control of breathing (part II of the model *via* term 1 of Equation 6: **V̇ ^{0}** [1 − C

^{REM}/2C

^{50}]). 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 Pco

^{2}increases. 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 ·[P

^{E}− P

^{E_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 Pco

^{2}. 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.

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̇− Pco^{2}response 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̇− Pco^{2}response 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½k^{e0}that 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½k^{e0}value 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½k^{e0}in 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 C^{50}(now only term 1 of Equation 6is operative). This occurs in a situation in which Pco^{2}remains 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 C^{50}(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 C^{50}value (0.7 ng/ml) for ventilation was measured at a fixed end-tidal Pco^{2}of 55 Torr.10The observed C^{50}is 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 C^{REM}> 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̇− Pco^{2}response 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 Spo^{2}, as endpoint. Low Spo^{2}values (< 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 Spo^{2}values 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 (Spo^{2}< 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 Spo^{2}may 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 Pco^{2}, ventilation at a fixed Pco^{2}of 55 Torr, and the V̇− Pco^{2}response 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 E^{max}function 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-E^{max}function: 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 E^{max}nature of term 1 and the multiplicative nature of the model, no apnea may be described or predicted. At two times the C^{50}, the Bouillon *et al.* 12model predicts ventilation levels more than 2 l/min, whereas in our approach, 2C^{50}equals the concentration at which apnea occurs (C^{APNEA}= 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 k^{e0}). 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 Pco^{2}yielded 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 C^{REM}, which yielded changes toward the steady state in both Pco^{2}and 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.

## References

^{2}and O

^{2}and 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

^{2}sensitivity. PLoS Med 2007; 4:e239

*versus*morphine on hypoxic and hypercapnic breathing in healthy volunteers. Anesthesiology 2003; 99:788–98