A biomarker of opioid-induced respiratory toxicity in experimental studies

Summary Opioids are commonly used painkillers and drugs of abuse and have serious toxic effects including potentially lethal respiratory depression. It remains unknown which respiratory parameter is the most sensitive biomarker of opioid-induced respiratory depression (OIRD). To evaluate this issue, we studied 24 volunteers and measured resting ventilation, resting end-tidal PCO2 (PETCO2) and the hypercapnic ventilatory response (HCVR) before and at 1-h intervals following intake of the opioid tapentadol. Pharmacokinetic/pharmacodynamic analyses that included CO2 kinetics were applied to model the responses with focus on resting variables obtained without added CO2, HCVR slope and ventilation at an extrapolated PETCO2 of 55 mmHg (V˙E55). The HCVR, particularly V˙E55 followed by slope, was most sensitive in terms of potency; resting variables were least sensitive and responded slower to the opioid. Using V˙E55 as biomarker in quantitative studies on OIRD allows standardized comparison among opioids in the assessment of their safety.

Opioids are drugs of abuse and produce potentially lethal respiratory depression (RD) Isohypercapnic ventilation is the most sensitive biomarker of opioidinduced RD Baseline variables are least sensitive biomarkers of opioid respiratory depression Comparative and experimental opioid studies should use isohypercapnic ventilation

INTRODUCTION
Opioids are commonly used in clinical practice to manage acute and chronic pain and suppress the sympathetic stress response during surgery. 1 However, opioids are toxic and possess several adverse effects including abuse, addiction and potentially life-threatening respiratory depression (RD). The combination of addiction and RD is the root cause of the current opioid crisis and numerous deaths in the US and other Western nations. [2][3][4] Opioids are effective and possibly even the most effective pain killers available. 1 Consequently and despite their known toxicity, their clinical use will continue and new opioids are being developed with the promise of reduced respiratory toxicity. 5,6 To validate whether opioids have an improved respiratory safety profile compared to other drugs of the same class, a comprehensive understanding of their respiratory effects on the human ventilatory control system is necessary. A large number of studies explored the effects of opioids on breathing, 1 but it remains unclear which of the many respiratory variables is the most sensitive and consequently most useful biomarker of the impact that opioids have on ventilatory control. In the current study, we address this issue by performing a physiology based pharmacokinetic/pharmacodynamic (PK/PD) analysis of the effect of an opioid on the hypercapnic ventilatory response (HCVR), integrating ''resting'' data, obtained before CO 2 inhalation (Figure 1), in the analysis. The model incorporates the CO 2 alveolar mass balance, [7][8][9] and describes the horizontal part of the curve where ventilation is independent of CO 2 . 10 The resultant is a function of ventilation versus steady-state PCO 2 . Four biomarkers were extracted from the analysis: resting ventilation ( _ V E ), resting end-tidal PCO 2 (P ET CO 2 ), the ventilatory recruitment threshold (VRT or the end-tidal PCO 2 at which ventilation begins to rise on CO 2 inhalation) and the slope of the HCVR. Furthermore, we calculated ventilation at an extrapolated end-tidal PCO 2 of 55 mmHg ( _ V E 55). 11,12 These five biomarkers were subsequently modeled over time using a sigmoid E MAX model (with parameters potency, C 50 , and hysteresis half-life, h½) to determine which one was most sensitive to the opioid. Experiments were performed in healthy and opioid-naïve male and female volunteers, who either received a low dose or a high dose of the opioid tapentadol.

RESULTS
Individual plasma concentrations and population data fits following low-and high-dose tapentadol ingestion are given in Figure 2. The PK data were analyzed with a two-compartment model, with the following parameter estimates: Volume of compartment 1 (V 1 ) = 90 G 6 (median G standard error of the estimate) L/70 kg, volume of compartment 2 (V 2 ) = 557 G 24 L/70 kg, elimination clearance from compartment 1 (CL) = 197 L/h (for a 70 kg individual) with between-subject variance 0.06 G 0.02, intercompartmental clearance (CL 2 ) 726 G 34 L/h, tapentadol absorption time from intake to blood 1.2 G 0.1 h and residual error 0.015 G 0.002.
Examples from one subject of the pharmacodynamic analyses are given in Figure 3 (the symbols are the measured data, the lines the estimated model output). It shows the analytical steps from inspired PCO 2 to the HCVR with its parameters resting _ V E , resting P ET CO 2 , VRT and the slope of the HCVR. Note that the HCVR is plotted versus the steady-state CO 2 partial pressure. In our model the steady-state PCO 2 is equivalent to the effect-site PCO 2 , i.e., the site where CO 2 is presumably sensed (see below). The extracted parameter estimates and calculated _ V E 55 were next analyzed using sigmoid E MAX models. The individual parameter estimates over time, the median values of the posterior distribution from the NONMEM Bayesian step G2.5 and 97.5 percentiles and estimated population fits are given in Figure 4 for the high-dose opioid dataset (panels A-E; all data presented are derived from the modeling analysis). These pharmacodynamic parameter estimates are given in Table 1. For all parameters, shape factor g was not different from 1 and therefore fixed to 1. Resting _ V E , resting P ET CO 2 and VRT had similar dynamics with similar values for potency (C 50 = 0.81 ng/mL, meaning that at this concentration their value drops by 50%) and hysteresis (h½ = 1.3 h), whereas the dynamics of parameters S and _ V E 55 were distinct from the other parameters (p < 0.01; Table 1). The picture that emerges from the analyses is that resting _ V E , resting P ET CO 2 and VRT respond slower with lesser sensitivity to the opioid than S and _ V E 55. Comparing S and _ V E 55, _ V E 55 was most sensitive in terms of potency with a C 50 value of 0.08 ng/mL, i.e., 10-fold more sensitive than resting parameters, versus C 50 for S of 0.23 ng/mL. S displayed a lesser hysteresis (h½ = 0.6 versus 0.9 h). The differences in dynamics between the HCVR, particularly S and resting parameters is exemplified in Figure 5, showing the effect of 200 mg tapentadol at t = 0 (control) and 2 and 4 h after opioid intake. The decrease in slope in this specific subject peaked at 2 h (and afterward S returned toward pre-drug values), whereas resting parameters had their peak at t = 4 h.
When comparing the opioid doses, the dose-response effect was largest for parameter S, whereas the smallest dose effect was observed for resting _ V E (in comparison, dose effect S > _ V E 55> VRT > resting P ET CO 2 > resting _ V E ; Figure 4F). This indicates that parameters derived from the HCVR were more sensitive to an increasing opioid dose compared to resting data.

DISCUSSION
It is important to evaluate and compare the respiratory effects of existing opioids and those in development to ensure their safe clinical use. There are different biomarkers of the respiratory effect of opioids, but their sensitivity and utility remain understudied. To quantify the effect of opioids on ventilatory control, some studies analyzed the effect of opioids on resting respiratory variables (e.g. respiratory rate) and estimated rather high C 50 values. This implies that the tested opioids had a low sensitivity or potency for RD and consequently are erroneously considered relatively safe. For example, Mildh et al. 13 studied the effect of two potent intravenous opioids, alfentanil and fentanyl, in eight healthy volunteers, and measured plasma concentration, resting ventilation, arterial PCO 2 and respiratory rate under poikilocapnic conditions, i.e., without strictly controlled end-tidal CO 2 . PK/PD modeling was used to estimate the opioid The response was obtained by performing steps in end-tidal PCO 2 . Each circle is a 1-min ventilation average. The orange cloud includes data obtained without any added inspired CO 2 (resting or poikilocapnic data). The blue cloud includes data obtained at inspired CO 2 concentrations and is the HCVR. iScience Article concentration that produces 50% decrease in resting _ V E , respiratory rate and increase in PCO 2 . The C 50values were 195 ng/mL for alfentanil and 3.5 ng/mL for fentanyl, with identical values for the three distinct biomarkers. We earlier estimated C 50 -values obtained for isohypercapnic _ V E (i.e., _ V E measured at a clamped end-tidal PCO 2 of 50-55 mmHg) in the order of 1.0 ng/mL for fentanyl and 50 ng/mL for alfentanil. 14,15 This indicates a 3-fold lower opioid sensitivity in data obtained at poikilocapnia compared to data obtained at isohypercapnia. The reason for the differences is that at poikilocapnia the opioid-induced reduction in _ V E and subsequent rise in arterial PCO 2 interact within the CO 2 chemoreflex feedback loop in such a way that the net changes in either parameter are minimized. 16 One might argue that the data obtained under resting or poikilocapnic conditions are more realistic, however, we would counter-argue that the model estimates depend on many experimental conditions that require consideration, such as the dose and the speed at which the opioid is administered. We showed that an opioid that is administered rapidly, causing a high peak in plasma and brain concentration, rapidly silences respiratory rhythm generation, whereas the same drug that is administered slowly allows accumulation of CO 2 that will offset (part of) the RD. 8 Modeling such effects without considering CO 2 kinetics will then result in a high C 50 value that holds the assumption of relative safety. Hence, we consider the C 50 -estimates derived from the opioid effect on isohypercapnic ventilation a more realistic approximation of the opioid's respiratory potency, i.e., not underestimating the opioid's potency. The current results provide additional arguments for this opinion.
We compared the pharmacodynamic properties of five biomarkers of ventilation, of which two were obtained before any CO 2 inhalation, resting _ V E and P ET CO 2 , and three indices were part of the HCVR curve, VRT, S and _ V E 55. Resting data together with VRT behaved similarly in terms of the opioid pharmacodynamics with similar values for C 50 and h½. S and _ V E 55 had distinct pharmacodynamics with some differences in their parameter estimates. Differences between the two sets of parameters were that resting data and VRT displayed a lesser potency and had a slower onset/offset in opioid effect than S and _ V E 55. The potency ratios were such that S and _ V E 55 were 3.5 and 10 times more potent than the other parameter set, whereas h½ differed by a factor of 2 with a more rapid onset/offset for S and _ V E 55. We relate the slower response and lesser opioid sensitivity of the two resting parameters and VRT to the interactive chemoreflex-related effect of arterial CO 2 on _ V E and vice versa, as discussed above. 16 The higher sensitivity of _ V E 55 relative to S, may be related to the fact that _ V E 55 is dependent on the slope as well as on the position of the HCVR curve ( Figure 5). The high opioid sensitivity observed for _ V E 55 indicates that this specific index is the most sensitive and, as we argue, useful biomarker when quantifying the impact of opioids on ventilatory control. The choice of ''best'' clinical opioid when comparing these drugs for safety during pain relief, however, is determined by their utility function or therapeutic ratio. 14, 17 In other words, opioids that have the largest separation between C 50 for analgesia and RD, under the condition that C 50 for analgesia < C 50 for RD, might be considered the opioid with theoretically a greater likelihood for analgesia than RD, although in the clinical setting this still depends on dose and speed of drug administration. We here define C 50 for analgesia as a 50% increase in ability to fend nociceptive stimuli (e.g. a 50% increase in pain pressure threshold) and C 50 for RD as a 50% depression of isohypercapnic ventilation or _ V E 55. Unfortunately, most opioids show greater or equal potency for RD than analgesia. For example, for oxycodone, C 50 s are 0.08 ng/mL for analgesia and 0.05 ng/mL for RD (unpublished observation), for fentanyl equivalent values are 1.8 and 1.0 ng/mL, 14 whereas for morphine C 50 was equal for both endpoints (10 ng/mL). 18 In all cases, the analysis of RD was based on isohypercapnic breathing (P ET CO 2 fixed at 50-55 mmHg) and the analysis of analgesia was based on experimental studies on antinociception. Hence, one cannot blindly extrapolate iScience Article these data to the clinical setting. Our analysis is particularly useful in experimental comparative studies of opioids and other drugs that affect ventilatory control, though.
We administered tapentadol as example opioid. In contrast to classical opioids, it has a dual mechanism of action. 19 It activates the m-opioid receptor and simultaneously inhibits neuronal reuptake of noradrenaline. This later mechanism is analgesic as noradrenaline inhibits post-synaptic nociceptive neurons in the spinal cord that express a 2 -adrenergic receptors. Because activation of the adrenergic system has little respiratory effects, we expect that tapentadol has a respiratory advantage over other opioids. However, the noradrenergic effect takes time to develop, very similar to the time it takes for selective antidepressants with the same mechanism of action to produce effective mood improvement (selective serotonin and noradrenaline reuptake inhibitors such as venlafaxine and duloxetine). Further studies in patients taking tapentadol for at least 2 weeks are needed to determine whether the separation of C 50 for analgesia and RD exceed that of, for example oxycodone. Traditionally, the HCVR is analyzed by fitting the ventilation data to the following formula _ V E = S 3 (P ET CO 2 -B), 16 where B is the extrapolated P ET CO 2 at zero ventilation (apneic threshold) and S equals per definition the combined sensitivity of the peripheral and central chemoreceptors. We plotted the apneic thresholds of a control HCVR curve and a curve at peak RD, as based on _ V E 55, in Figure 5. When considering the magnitude of the opioid effect, the leftward shift of the apneic thresholds gives little information without any knowledge on the change in slope. Biomarker _ V E 55 captures both and hence is our choice rather than just S or B for quantifying drug effect on ventilatory control.
The HCVR is a function of PCO 2 (Figure 1). In steady-state hypercapnic experiments, in which inspired CO 2 is raised and P ET CO 2 and _ V E are at steady state, i.e., both variables show no further increase over time, we assume that P ET CO 2 closely approaches brain tissue or effect-site PCO 2 . Under these circumstances, the HCVR is measured at the site of the central chemoreceptors, assuming that the ventilatory response to CO 2 arises predominantly from the activation of central chemoreceptors. We do stress that the effectsite PCO 2 does not correspond to the anatomical sites of CO 2 sensing. The sites of central CO 2 sensing are still under debate and the mechanism of CO 2 sensing depends on complex mechanisms not captured in our pharmacodynamic model. Steady-state conditions are not met when non-steady-state experiments are conducted or experiments in which CO 2 rises without any steady state in P ET CO 2 or _ V E . 20 The resultant HCVR is then not a function of the effect-site or steady-state PCO 2 . However, if the analysis considers CO 2 kinetics, the response will be a function of the effect-site or steady-state PCO 2 irrespective of the method of inducing hypercapnia. Our analysis is therefore a major advantage over other methods as it does not require an experimental steady state in _ V E or P ET CO 2 .
We modeled _ V E as linear function of PCO 2 beyond the VRT (Figures 3 and 5). Others used a non-linear exponential function to describe the _ V E -PCO 2 relationship and included resting data. 7,9 It is our consistent iScience Article observation that the HCVR over the PCO 2 range measured (beyond the VRT) is best approximated by a linear function (Figures 1, 3, and 5). Moreover, the dynamics of the initial horizontal part of the curve is dissimilar compared to the dynamics of the rest of the curve, as we demonstrated here, and therefore requires a distinct approach. We are aware that at higher PCO 2 values the curve flattens and the complete HCVR has a sigmoid shape. We argue that our approach best reflects the physiological changes that occur following inhalation of up to 9% CO 2 . 21

Conclusions
We systematically analyzed the effect of a potent opioid on five biomarkers of the HCVR and observed that _ V E 55, calculated from the HCVR curve, is the most sensitive quantitative index of opioid-induced RD. In addition, we showed that resting variables such as resting _ V E and P ET CO 2 respond slower and with a lesser sensitivity to the opioid than the slope of the HCVR and _ V E 55 and hence underestimate the opioid respiratory effect. We argue that _ V E 55 is the appropriate index to be used in experimental studies on the effect of opioids (and other drugs) on ventilatory control. Adoption of this index as biomarker of opioid effect on ventilation will reduce the rather diverse outcome of opioid respiratory studies and will allow a clear quantitative comparison among opioids in the assessment of their safety.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:   Differential effect of resting data and data obtained at an elevated PCO 2 observed in subject 001 and derived from the data modeling. The black line is the HCVR curve obtain before any opioid intake (time 0). The horizontal part of the curve is the so-called dog-leg and reflects resting data without added CO 2 (poikilocapnic data). The line deflects upwards at the ventilatory recruitment threshold and turns into the response obtained at added inspired CO 2 (isohypercapnic data). The next response, red line, is obtained 2 h after 200 mg oral tapentadol intake. The slope of the HCVR has decreased to its lowest value measured, whereas the resting ventilation has not changed from control, indicative of the greater sensitivity of the response obtained at added inspired CO 2 than resting data obtained at poikilocapnia. At t = 4 h after opioid intake, green line, the slope returns toward control, value, whereas resting ventilation now shows its peak depression effect, indicative of the slower response of resting data compared to data obtained at isohypercapnia. On the xaxis, the socalled apneic thresholds (or PCO 2 at apnea) for the control curve and the one obtained at t = 2 h are depicted as black and red symbols, respectively. The rising part of the curve is commonly analyzed by the function Ventilation = S 3 (end-tidal PCO 2 -B) where B is the apneic threshold.

Materials availability
This study did not generate new unique reagents.
Data and code availability d All data reported in this paper will be shared by the lead contact upon reasonable request.
d The NONMEM code reported in this paper will be shared by the lead contact upon reasonable request. d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon reasonable request.

Study conduct and ethics
The current report is part of a larger project in which two opioids (oxycodone and tapentadol) were tested; here we present the results of the study on tapentadol. The study was performed from January 2019 to December 2020 after approval of the study was obtained from the Medical Research Review Board Leiden-the Hague-Delft (METC-LDD, located in Leiden, the Netherlands) and the Central Committee on Research Involving Human Subjects (the Hague, the Netherlands) on November 28, 2019. All subjects gave written informed consent prior to enrollment in the trial.

Study design and population
We performed a single-center, randomized study on the effect of an opioid on ventilation as measured by the steady-state hypercapnic ventilatory response (HCVR). We studied 24 adult healthy individuals (12 men/12 women; mean age 23 years with range 19-33 years, body mass index 23.4 kg/m 2 with range 19.1-28.7 kg/ m 2 , and weight 74 kg with range 57-103 kg). A 20 G arterial line was inserted in the radial artery of the left or right arm for blood sampling. Subjects received a low (100 mg) or high oral dose (200 mg) of the opioid tapentadol immediate release (Grü nenthal GmbH, Aachen, Germany). Each subject only received one dose, and the order of study drug allocation was random and double blind. The drug was ingested with 100 mL tap water after a control HCVR was obtained. Thereafter HCVRs were obtained at 1-h intervals for 8 hours.

Study endpoints/outcome measures
To measure ventilation ( _ V E ), the participants inhaled a gas mixture via a face mask coupled to a pneumotachograph and pressure transducer system (Hans Rudolph Inc., USA). iScience Article CO 2 and N 2 delivered by three mass flow controllers (Bronkhorst BV, the Netherlands) that were controlled by custom-made software (Leiden University). The software allowed the collection of ventilation data on a breath-to-breath basis and enabled variations in inspired gas concentrations to achieve the desired endtidal oxygen and carbon dioxide concentrations (P ET O 2 and P ET CO 2 ), independent of the ventilatory response. The concept of this dynamic end-tidal forcing approach is described elsewhere in detail. 22,23 In brief, in this study, subjects initially inhaled a normoxic gas mixture (P ET O 2 110 mmHg) without any added inspired CO 2 for 5-6 min, after which we applied three to four steps in P ET CO 2 with a fixed sequence: 4.5, 9, 13.5 and 18 mmHg above resting P ET CO 2. In case the procedure tended to last longer than 20 min for the first three CO 2 steps, the last step (18 mmHg) was omitted. The duration of each step was 6-8 min and ended when we observed a steady state in _ V E (at least 2 minutes no further increase in _ V E ). In-between HCVRs the subjects breathed room air. Inspired and expired O 2 and CO 2 were measured at the mouth (Datex Capnomac, Finland), and for monitoring purposes the ECG and oxygen saturation (Datex Cardiocap, Finland) were measured throughout the experiment day.
Blood samples were drawn from the arterial line for measurement of the tapentadol plasma concentrations. A 4 mL blood sample was drawn before dosing and at 1-h intervals for 7-h. The samples were analyzed by Ardena Bioanalytical Laboratory (Netherlands). Tapentadol concentrations were measured by tandem mass spectrometry (LC-MS/MS) methods, validated over a range of 0.1 to 1,900 ng/mL, in one batch. Within-run precision (% coefficient of variation) and accuracy (% bias) was 4.0 and 0.7 at the lowest level of quantitation, respectively, and 0.6 and 2.8 and the highest level, respectively.

Sample collection
A formal power analysis was not performed but we enrolled a convenience sample of 24 subjects with 12 subjects per dose, which is common practice in similar PK/PD studies. 24 Subjects were randomized to a low or a high dose regimen using computer-coded randomization list by an independent statistician. which was made available to the pharmacy who prepared the drug in unmarked syringes.

QUANTIFICATION AND STATISTICAL ANALYSIS
We performed a physiology-based pharmacokinetic/pharmacodynamic (PK/PD) analysis to determine the effect of the opioid on the different components of the HCVR: resting ventilation, resting P ET CO 2 , slope of the HCVR and _ V E at an extrapolated P ET CO 2 of 55 mmHg ( _ V E 55). The pharmacokinetics and pharmacodynamics of were analyzed with NONMEM VII (Icon Plc., USA), a software package for nonlinear mixed effects modeling, using a population approach. 25 The pharmacokinetic data were analyzed first, using a two-compartmental model. Parameter estimation was done using Stochastic Approximation Expectation Maximization (SAEM), followed by objective function evaluation using Importance sampling (IMP). Rather than using an absorption compartment, we used the zero-order infusion duration from tablet into the blood. The adequacy of the pharmacokinetic model was evaluated by inspection of the individual data fits and the goodness-of-fit plots: measured concentrations versus individual predicted concentrations, measured concentrations versus population predicted concentrations, conditional weighted residuals versus time, and normalized prediction discrepancy error versus time. No covariates were included in the pharmacokinetic model analysis. The output of the PK analysis was used as one of the input variables to the pharmacodynamic analysis.
The _ V E data were averaged over 1-min and these 1-min averages were used in the data analysis. The complete _ V E dataset obtained during one HCVR run was analyzed to obtain the HCVR curve ( Figure 1). To that end, we defined various components of the response that were included in the analysis: resting ventilation ( _ V B ), resting P ET CO 2 , the ventilatory recruitment threshold (VRT) or the P ET CO 2 at which ventilation starts to rise linearly, and the slope of the HCVR (S). Resting data were obtained prior to any CO 2 inhalation. We also calculated for each HCVR curve the extrapolated _ V E at a P ET CO 2 of 55 mmHg ( _ V E 55), as an additional measure of drug effect. Finally, we quantified the hysteresis between plasma opioid concentration and effect (blood-brain equilibration half-life or h½) and opioid potency (C 50 or the opioid concentration in the brain causing a 50% effect); h½ and C 50 were determined for all parameters, _ V B , resting P ET CO 2 , VRT, S and _ V E 55 independently. By including the hysteresis in our models, the HCVR becomes a function of the effect-site PCO 2 or PCO 2 at the site were CO 2 is sensed.

OPEN ACCESS
iScience 26, 106520, April 21, 2023 11 iScience Article The relationship between CO 2 content (C) and its partial pressure (P) was assumed to be linear, so that P = l 0 3 C, where l 0 = 0.115 mmHg/(mL CO 2 in 100 mL blood). The following CO 2 alveolar mass balance equations were used for the lungs and body (approximating the body by one compartment): 8,21 V ALV $ d À arterial P CO 2 Á dt = ð _ V E À _ V D Þ $ inspired P CO 2 À arterial P CO 2 + l 1 $ _ Q$ À venous P CO 2 À arterial P CO 2 Á and V TS $ d À venous P CO 2 Á dt = _ Q $ À venous P CO 2 À arterial P CO 2 Á + l 2 $ _ V CO 2 ; Where V ALV is the alveolar volume, the arterial carbon dioxide pressure is assumed to equal alveolar pressure, _ V D dead space ventilation, _ Q cardiac output, vP CO2 the venous carbon dioxide partial pressure, V TS the apparent CO 2 distribution volume in tissue, _ V CO2 the carbon dioxide production, l 1 = q 3 P BW /l 0 /100 z 10 and l 2 = 100 3 l 0 , where q = 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. _ V ALV was fixed to 3 L and _ V D to 1.8 L/min; V TS and _ Q were estimated to be 10 and 4 L, respectively.
_ V E was assumed to depend on brain tissue carbon dioxide pressure (brain tissue P CO2 ) as: t $ d À brain tissue P CO 2 Á dt = arterial P CO 2 À brain tissue P CO 2 _ V E = baseline _ V E + S $ H$ À brain tissue P CO 2 À VRT Á and HðxÞ = d $log where t is a time constant (estimated to be 2.5 min) 22,23,26 and H is the so-called ''hinge'' function (Columbia University, a continuous hinge function for statistical modeling: https://statmodeling.stat.columbia.edu/ 2017/05/19/continuous-hinge-function-bayesian-modeling/), with d fixed to 0.1 and x = (brain tissue PCO 2 -VRT). For each respiratory run, parameters resting _ V E , resting P ET CO2, VRT, S were estimated or calculated ( _ V E 55), which were then used to estimate h½ and C 50 using sigmoid E MAX models. The sigmoid E MAX model had the form: EffectðtÞ = ðEffect at baselineÞ3½Parameter valueðtÞ=Parameter valueðtÞ + C 50 g where g is a shape parameter, where Effect at baseline is the control response prior to drug intake. P-values <0.01 were considered significant.

ADDITIONAL RESOURCES
Registration of the study was at the Dutch Clinical Trial Register, available at https://clinicaltrialregister.nl/ en/trial/23681 (id 23681) on December 31, 2018.