#### What We Already Know about This Topic

#### What This Article Tells Us That Is New

^{−1}· min

^{−1}

^{1,2}As 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).

^{3}Especially, 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.

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

^{1}However, 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}

*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–13}In 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}

Fig. 1 Image Tools |
Equation 8 Image Tools |

^{9,10}Although 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.

^{14}Hence, 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 breathing

^{3}; (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).

^{10}In 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

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

*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.

^{17}We 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

_{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.

^{10}

##### Data Analysis

##### Remifentanil Pharmacokinetics.

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

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

_{0}× C, where λ

_{0}= 0.863 Torr/(ml CO

_{2}in 100 ml blood) (fig. 2).

^{19}The following mass balance equations were used for the lungs and body (approximating the body by one compartment):

Equation (Uncited) Image Tools |
Equation 1 Image Tools |

_{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̇CO2**was estimated from the baseline end-tidal carbon dioxide pressure and V̇. These baseline values,

**Q̇**and V

_{TS}were parameters to be estimated.

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

_{2}(P

_{E}) on ventilation under hyperoxic conditions can be modeled as follows (fig. 2)

^{14–16}:

*et al.*

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

_{0}is the apneic threshold when C

_{REM}= 0 and C

_{100}the concentration remifentanil that causes a doubling of B. Rewriting Equation 3 and 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:

**V̇0**) = G · (P

_{E_0}− B

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

_{50}=

Equation 5 Image Tools |
Equation (Uncited) Image Tools |

_{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.

_{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:

_{E}is the measured value and

**P̂E**the predicted value. The variance of ε should be smaller (ς

_{A}

^{2}) when P

_{E}>

**P̂E**and larger when (ς

_{B}

^{2}) when P

_{E}<

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

**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.

_{2}measurements are inaccurate.

##### Modeling Minute Ventilation.

_{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.

^{18}The 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.*

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

^{20}The 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

^{−1}· min

^{−1}. The infusion was terminated when the effect site concentration had reached 5 ng/ml.

#### Results

^{−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

*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 4 are in agreement with this statement. The pharmacokinetic parameter estimates were not significantly different from those of Minto

*et al.*

^{17}except 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

^{−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

_{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

Fig. 7 Image Tools |
Fig. 8 Image Tools |

^{−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 8 for 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.

_{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

*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,12}Although 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

_{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,2}Part 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 6 is 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.

_{2}response slope.

^{10,21,22}However, 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.

^{23}As 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}

_{e0}that we observed (0.5 min) is in the same range as values from Babenco

*et al.*

^{8}using 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–26}However, this is only true when the drug is not overdosed. Otherwise, severe respiratory depression with apnea should always be in mind.

^{4}Our data further indicate that the chances of apnea are increased when the subject is asleep (BIS values = 80) with propofol.

_{50}(now only term 1 of Equation 6 is 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.*

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

^{10}The observed C

_{50}is a low value, probably related to a differences in parameterization (see Eq. 4 of Ref.

^{10}) and the fact that we analyzed the whole remifentanil–propofol surface rather than just the remifentanil-ventilation relationship.

_{REM}> 0 = 1.6 l · min

^{−1}· Torr

^{−1}; fig. 1).

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

^{15}Furthermore, 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.

^{14}The value observed by us is in agreement with the value presented by Babenco

*et al.*

^{8}after 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}

##### Apnea

*et al.*

^{27}comes 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.*,

^{27}we 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,16}This 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

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

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

^{3,12}We 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

*et al.*

^{29}modeled 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

*et al.*

^{11,12}were 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.*

*et al.*

^{11–13}used 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}):

_{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.*

^{12}model 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.*

^{12}to 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.*

^{12}suggest 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).