FormalPara Key Points

Weak pharmacokinetic (PK) interactions are expected when the normal doses of opioids are concomitantly taken with overdosed benzodiazepines.

Pharmacodynamic interactions may play a more important role than PK interactions in causing drug–drug interactions between opioids and benzodiazepines.

1 Introduction

A key finding in clinical pharmacology and therapeutics is that most overdose fatalities involve multiple drug classes, complicating the safety of a specific drug. In fact, the coadministration of drugs is likely to be a risk factor. For example, sedatives were estimated to be involved in 11,843 and 1847 deaths in 2014 and 1999, respectively, whereas sedatives were virtually never solely implicated in those deaths [1]. Opioids are drugs that can act on opioid receptors and produce morphine-like effects. They have been widely used for pain relief for many years. However, opioids have both benefits and overdose side effects, including nausea, vomiting and coma [2, 3]. In the past two decades, opioid prescribing has increased tremendously in the USA. In total, 16,651 deaths were related to opioid medications in 2010 [4]. For example, the percentage of poisoning deaths caused by heroin itself increased 42.6% between 2007 and 2014 [1]. Opioids are also one of the most frequently used drugs in drug combinations. Despite the therapeutic benefits of opioids, their safety issues cannot be ignored. The combination of opioids and other drugs can cause even greater risks. For example, the percentage of poisoning deaths increased to 97.2% when heroine was combined with other drugs during 2007–2014 [1].

Benzodiazepines are one of the most commonly coadministered drugs and are often prescribed for patients with anxiety disorders, muscle spasms or major depression [5]. Since the 1970s, many researchers and physicians have investigated the coadministration of opioids and benzodiazepines [6]; approximately 5000 publications related to these two classes were published between 1970 and 2012 [6]. Studies have indicated that, although the risk of overdosing on benzodiazepines alone was mild, the combination of opioids and benzodiazepines (especially overdosed benzodiazepines) posed a potential danger to patients because of the increased risk of synergistic respiratory depression and overdose death [3, 7,8,9,10]. Research has confirmed that drug–drug interactions (DDIs) between opioids and benzodiazepines do exist, and many studies have attempted to illustrate how these two drug classes interact with each other [6, 11, 12].

Two basic factors can contribute to DDIs: pharmacokinetic (PK) and pharmacodynamic (PD) interactions, with PK interactions the most common. These may occur when a coadministered drug causes a change in the absorption, distribution, metabolism and/or excretion (ADME) of another drug [13]. Another mechanism that underlies drug interactions can be explained from the PD perspective. This kind of interaction occurs when two or more drugs are involved in the mechanisms of action that result in the additive, synergistic or opposite physiological outcome.

It is believed that benzodiazepines can alter the PKs of opioids [6, 14]. Opioids primarily undergo phase I metabolism through cytochrome P450 (CYP) 3A4 enzymes and, therefore, may have significant interactions with other coadministered drugs that are also sCYP3A4 substrates, inhibitors or inducers [15]. Meanwhile, some benzodiazepines have been reported as CYP3A4 inhibitors since they are also mainly metabolized by the CYP3A4 enzyme [11, 14, 16,17,18,19]. Some studies suggested that coadministration of benzodiazepines with opioids can potentially increase opioid exposure. For example, research utilizing human liver microsomes demonstrated that midazolam was a moderate mechanism-based inactivator of buprenorphine N-dealkylation, which could cause time- and concentration-dependent inhibition of norbuprenorphine, a metabolite of buprenorphine produced mainly by the CYP3A4 enzyme [18]. By quantitatively analyzing the plasma concentrations of oxycodone and clonazepam, Burrows et al. [20] reported that concomitant clonazepam intake could reduce the metabolism of oxycodone. However, almost all these study cases were performed using rat or human liver microsomes. Therefore, we ask whether benzodiazepines can affect the metabolism of opioids in the human body. In this work, we used physiologically based PK (PBPK) modeling to address this issue.

To study the DDIs between opioids and benzodiazepines, we constructed both full and minimal PBPK models for a set of commonly used opioids and benzodiazepines. Three opioid drugs (fentanyl, oxycodone and buprenorphine) and four benzodiazepine drugs (alprazolam, diazepam, midazolam and triazolam) were selected for the PK interaction studies. We used molecular docking scores to calculate the missing inhibition parameters, Ki, and conducted sensitivity analysis to evaluate the impact of Ki on the PBPK modeling results. A PBPK model is useful only after it is validated with observed data, such as the drug’s concentration–time (CT) curves. The DDI profile of a drug combination is then generated through the simulation.

2 Methods

2.1 Representative Opioids and Benzodiazepines

We selected the opioids oxycodone, buprenorphine and fentanyl to study the PK interactions. CYP3A4 enzyme is an important metabolizer of these three drugs [21]. Oxycodone is mainly metabolized by CYP3A4, and noroxycodone is its metabolite, although a small fraction of oxycodone also goes through CYP2D6 metabolism to produce oxymorphone [22]. Buprenorphine is mainly metabolized by CYP3A4, CYP3A5 and CYP2C8 [23, 24]. The phase II metabolism by uridine diphosphate glucuronosyltransferase (UGT) is almost negligible for oxycodone and fentanyl since it contributes very little to their metabolism [25], although the metabolism of buprenorphine by UGT1A1, UGT1A3, UGT1A4 and UGT2B7 was taken into account [26] in this study.

We selected four frequently used benzodiazepines (alprazolam, diazepam, midazolam and triazolam). Alprazolam, midazolam and triazolam are predominately metabolized by CYP3A4 and CYP3A5 enzymes [27, 28]. Specifically, diazepam is mainly metabolized by CYP3A4 and CYP2C19 to its metabolites, temazepam and nordazepam [29]. CYPC18 and CYP2C9 enzymes also make very limited contributions to its metabolism. Generally, all four benzodiazepines primarily undergo phase I metabolism by CYP enzymes, and phase II metabolism is negligible.

2.2 Physiologically Based Pharmacokinetic (PBPK) Modeling

We used SimCYP Simulator (v. 17, release 1; Sheffield, UK) software to build the PBPK models by conducting virtual clinical trials using 100 virtual healthy volunteers (“Sim-healthy volunteers”). Full PBPK models with multiple compartments representing different organs of the human body were applied to the model development for diazepam and all opioids. Meanwhile, minimal PBPK models with the consideration of only three or four compartments (with or without the single adjustable compartment) were developed for all benzodiazepines except diazepam. The software’s advanced dissolution, absorption and metabolism (ADAM) model was used to simulate the absorption process, particularly for diazepam. The ADAM model considers the complicated process of drug absorption and the interplay with the underlying physiological characteristics of the gastrointestinal tract [30, 31]. Figure 1 shows the schematic diagrams of a generic full PBPK model and an ADAM model, and Fig. 2 shows that of a minimal PBPK model. PK data for ADME from in vitro or in vivo experiments were collected to generate the “compound profiles” of the drugs in study. The inhibitor profiles of alprazolam, midazolam and triazolam were created based on their corresponding substrate profiles stored in the SimCYP drug library.

Fig. 1
figure 1

The PBPK model (left panel) and ADAM model (right panel). ADAM advanced dissolution, absorption and metabolism, PBPK physiologically based pharmacokinetic, PO oral

Full size image
Fig. 2
figure 2

Schematic structure of minimal PBPK model. CLin/CLout clearance into/out of single adjusting compartment, Fa fraction absorbed from gastrointestinal tract, Fg gut availability (fraction of drug escaping from the gut availability), GI gastrointestinal, Ka absorption rate, PBPK physiologically based pharmacokinetic, PO oral, QH blood flow rate from liver to systemic blood, Qpv blood flow rate from systemic blood to portal vein or from portal vein to liver

Full size image

2.2.1 PBPK Model Development and Verification

To conduct a PBPK simulation, both the system-related and drug-related parameters are needed. All the system-related parameters, including comprehensive physiological parameters such as blood rates and volumes of different organs, come from the SimCYP internal database. We performed PBPK simulations using a virtual population of 100 healthy adults aged 18–50 years. SimCYP applies an algorithm to incorporate known variabilities in those system-related parameters so that the virtual population is representative. The drug-dependent parameters come from literature or public databases (such as the PubChem database https://pubchem.ncbi.nlm.nih.gov/) when available, otherwise they were estimated by the SimCYP software. To validate the PBPK models, observed PK data, especially the CT curves, were collected to be compared with the results from simulations. The dosages of opioids and benzodiazepines were determined so that in silico modeling was consistent with the corresponding experimental settings. The drug-dependent parameters, including the PK parameters of both opioids and benzodiazepines, are listed in Table S1.

2.2.2 Inhibition Parameters

As per the previous research, we hypothesized that benzodiazepines were competitive inhibitors of opioids. When no experimental value of the inhibitory constant, Ki, of a metabolic reaction was available, we applied a molecular modeling approach to estimate its value. We first performed docking simulation using the Glide program implemented in Schrodinger’s small-molecule drug discovery suite (http://www.schrodinger.com) and then calculated the Ki value using Eqs. 2–3.

$$S\begin{array}{*{20}c} {k_{\text{f}} } \\ \leftrightarrow \\ {k_{\text{r}} } \\ \end{array} {\text{ES}}\mathop \leftrightarrow \limits^{{K_{\text{cat}} }} E + P,$$
(1)
$$\Delta G^{0} = - RT\ln K_{\text{eq}} = - RT\ln \frac{{k_{\text{f}} }}{{k_{\text{r}} }},$$
(2)
$$K_{i} = \frac{\left[ E \right]\left[ I \right]}{{\left[ {EI} \right]}} = \frac{{k_{\text{r}} }}{{k_{\text{f}} }},$$
(3)

where E, S, ES and P denote enzyme, substrate and enzyme-substrate complex, and product, respectively. kf is the forward reaction rate constant of E + S, and kr is the reverse reaction constant describing the rate of falling apart to E + S from ES. kcat is the forward rate constant of the formation of E + P. [E], [I] and [EI] are the concentration of enzyme, inhibitor and the complex. \(\Delta G^{0}\) is the binding free energy that can be estimated using the Glide docking score (kcal/mol). Keq = Kf/Kr, is the equilibrium constant for the reversible reaction. The calculated Ki value for each benzodiazepine is listed in Table S1.

2.2.3 Simulation of Opioid–Benzodiazepine Drug–Drug Interaction (DDI)

Following the validation of the PBPK models of opioids and benzodiazepines, we then conducted simulations of DDIs between two kinds of drugs using the built-in “Sim-healthy volunteers” virtual population. We retained the doses and formulations of three opioids—oxycodone, fentanyl and buprenorphine—to be consistent with the model validation processes but unified the dosing strategy of four benzodiazepines (alprazolam, diazepam, midazolam and triazolam). First, the single oral 10-mg dosage of four benzodiazepines was considered as a normal dose. This dosage was then changed to 500 mg (overdosed) and applied to these four drugs. Ten trials with ten subjects per trial, with an age range of 18–50 years, were simulated to evaluate the inhibitory effects of benzodiazepines on the exposure of opioids. Finally, the predicted DDIs of opioid–benzodiazepine drug combinations were compared for the normal dose and overdose scenarios.

2.2.4 Sensitivity Analysis for K i

Since the docking scores may not be accurate enough, we conducted a sensitivity analysis to explore the impact of Ki value on the DDI effects by changing the Ki value from one-tenth to 100-fold the current value. The sensitivity analysis was conducted using the built-in sensitivity analysis function in SimCYP.

3 Results

3.1 Validation of PBPK Models for Opioids and Benzodiazepines

The PBPK model of each drug was developed using in vitro and in vivo data from published literature or databases. Figure 3 shows the predicted and observed mean plasma CT profiles of all opioids and benzodiazepines. PK parameters, including area under the curve (AUC), maximal concentration (Cmax) and time to Cmax (tmax) are listed in Table 1.

Fig. 3
figure 3

Predicted concentration profiles of sublingual buprenorphine 4 mg, intravenous fentanyl 0.1 mg/kg, oral alprazolam 2 mg, oral midazolam 15 mg and oral triazolam 0.25 mg versus their observed data, respectively. Red open circle, blue open square and yellow open triangle represent the observed data. Black curves represent concentration–time curves, and gray dashed curves represent 95% confidence intervals of the population-based simulation of concentrations

Full size image
Table 1 Measurements for buprenorphine, fentanyl, alprazolam, midazolam and triazolam in 24 or 12 h
Full size table

The predicted AUC, Cmax and tmax of the two types of drugs were all within the standard deviation (SD) ranges of observed values, except for buprenorphine and fentanyl. For buprenorphine, the calculated AUC and tmax were slightly smaller than the lower bound values, albeit the calculated Cmax was very close to the observed one. No measured PK parameters were available for fentanyl. Encouragingly, as shown in Fig. 3, the observed CT values of these drugs were also within the confidence interval (CI) ranges (the upper and lower gray dashed lines) of the simulated population-based CT curves. Therefore, our PBPK models were well-validated.

3.2 DDI Prediction

The validated PBPK models of opioids and benzodiazepines were then applied to predict the DDIs between any opioid and benzodiazepine. The predicted AUC ratios and the Cmax ratios of DDIs for both scenarios, normal dose and overdosed benzodiazepines, are listed in Table 2. Figure 4 shows the results of the simulated DDI effect of different doses of benzodiazepines on three opioids, i.e., oxycodone, fentanyl and buprenorphine.

Table 2 AUC and Cmax ratio of DDI profiles in 24 h for normal dosage of opioids, including oxycodone 30 mg, buprenorphine 4 mg and fentanyl 0.1 mg/kg, and benzodiazepines, including alprazolam, diazepam, midazolam and triazolam with normal (10 mg) and overdose (500 mg), respectively
Full size table
Fig. 4
figure 4

The area under the curve ratio of oxycodone, buprenorphine and fentanyl with the presence of normal dose and overdose of four benzodiazepines in 24 h. AUC area under the curve, ND normal dose, OD overdose

Full size image

As shown in Table 2 and Fig. 4 for oxycodone and fentanyl, the fold changes of AUCs increased significantly for the scenario in which benzodiazepines were overdosed. Conversely, for buprenorphine, there was no obvious change in AUC ratios even when the benzodiazepines were coadministered at 500 mg. Note that the change in AUC ratios for fentanyl with normal or overdoses of benzodiazepines was significantly larger than that for oxycodone and buprenorphine, especially when alprazolam was considered. Even for the drug combination of fentanyl and alprazolam, the largest AUC ratio is only 1.98, suggesting PK interactions make no or limited contributions to DDIs between the two types of drugs. We also calculated the contributed portions of CYP3A4 to the total clearance (fmCYP3A4) of oxycodone, buprenorphine and fentanyl. The values of fmCYP3A4 for the three drugs were 34.49%, 1.47%, 92.63%, which better explains the DDI simulation results.

3.3 The Impact of K i Value

To estimate the impact of Ki value on the simulation results for each benzodiazepine, we conducted sensitivity analysis by adjusting the Ki value from the 1/10- to 100-folds of the predicted value. Specifically, we adjusted the Ki values of alprazolam, diazepam, midazolam and triazolam to be 0.163–163, 0.165–165, 0.217–217, 0.102–102 µM to investigate how AUC ratios of each opioid would change. The sensitive analysis results for the overdosed benzodiazepine scenario are summarized in Fig. 5. The predicted AUC ratios were consistent with the results shown in Fig. 4, and maximal AUC ratios were observed when applying 1/10 of the Ki values in DDI simulations. However, the impact of Ki values on DDI simulation result is limited, and the maximum AUC ratio is still no larger than 3.0.

Fig. 5
figure 5

The area under the curve ratio of oxycodone, buprenorphine and fentanyl with the presence of overdose benzodiazepines (500 mg) in 24 h when applying 1/10- to 100-fold Ki value to the drug–drug interaction models in 24 h. AUC area under the curve

Full size image

4 Discussion

The purpose of this study was to explore the PK interactions between opioids and benzodiazepines and to investigate the effects of benzodiazepines on the metabolism of opioids in the human body. Although previous studies have tried to uncover evidence of PK interactions between these two kinds of drugs, most were performed in in vitro settings and/or with animals. Because experimental DDI studies in human bodies are rare, our state-of-the-art simulation provides an alternative approach to study PK interactions between these two types of drugs.

The simulated concentration profiles of three opioids and four benzodiazepines were mostly well-predicted compared with the observed data, as shown in Fig. 3 and Table 1, except for the AUC ratio and tmax of buprenorphine, for which the predicted values were slightly out of the SD range of the observed data. Our simulation results suggested there were no PK interactions between opioids and normal doses of benzodiazepines. However, the PK interactions could be observed when concomitantly taking a normal dose of opioids and with severe overdoses of benzodiazepines. This observation is particularly true for fentanyl, which demonstrated the largest AUC ratio increase when overdosed benzodiazepines were administered concomitantly (Fig. 4 and Table 2). This finding may be because the predicted liver fmCYP3A4 value for fentanyl, 92.63%, was much higher than that for oxycodone (34.49%) and buprenorphine (1.47%). There were almost no PK interactions between buprenorphine and all four benzodiazepines, even when the benzodiazepine doses were set to 500 times the normal dose, because the fmCYP3A4 of buprenorphine was only 1.47% in the liver.

Among the four benzodiazepines studied, alprazolam had the strongest inhibitory effect on fentanyl (Fig. 4 and Table 2). The AUC ratio of fentanyl with coadministration of overdose alprazolam was 1.98, which is dramatically higher than the AUC ratios when fentanyl was coadministered with the other benzodiazepines. As such, the combination of fentanyl and overdosed alprazolam should be avoided, or close attention should be paid to adverse effects if the drugs are coadministered in clinical practice. To further validate this finding, we searched the literature and found that, during 2012–2017, one of the most frequent co-intoxicants reported with fentanyl was alprazolam (15.5%) [32].

In this study, we also applied molecular modeling techniques to solve the “missing parameter” problem in PK modeling and simulation. PK modeling is becoming an essential part of drug discovery and development, but it is very challenging to build predictive models since many PK parameters are not available in the literature or public databases. Using molecular modeling to assign parameters for PK modeling provides a practical way to investigate DDIs for drugs lacking experimental PK parameters. Furthermore, molecular modeling is not restricted to describe a drug’s metabolic process; it can also be used to study DDIs from the PD and poly-pharmacological perspectives.

This study has some potential limitations. First, the transporter DDI was not considered because of the lack of clinical evidence. However, some transporters in the blood–brain barrier (BBB) might cause significant interactions between the two types of drugs. Second, the inhibitory abilities of metabolites of benzodiazepines were ignored since we assumed the concentration of most of the metabolites was much lower than the parent drugs, leading to underestimated inhibitory effect of benzodiazepines on the opioids. Last, the molecular docking approach might not be accurate, and the predicted Ki value may differ from the real value. More advanced binding free energy calculation methods should be applied in future studies. Fortunately, our sensitivity analysis suggested that the predicted AUC ratios due to DDI were not sensitive to the Ki value. Therefore, the results of our DDI simulation should still be reliable.

5 Conclusion

The PK interactions between opioids (oxycodone, buprenorphine and fentanyl) and benzodiazepines (alprazolam, diazepam, midazolam and triazolam) are negligible when both drugs are coadministered with the normal doses. Weak PK interactions are expected when the normal doses of opioids are concomitantly taken with overdosed benzodiazepines. PD interactions may play a more important role than PK interaction in causing DDIs between opioids and benzodiazepines. However, variance still exists in terms of PK results when the three opioids are taken together with different benzodiazepines. In particular, obvious PK DDI occurs with concurrent use of fentanyl and overdosed alprazolam. Therefore, in daily clinical practice, alprazolam should be replaced with other benzodiazepines in patients addicted to fentanyl.