1. Introduction
1.1. GABAA Receptors
Gamma-aminobutyric acid type A receptors (GABAARs) represent a diverse family of proteins found throughout the nervous system and in non-neuronal tissues. These receptors assemble into transmembrane channels, either as homomers or heteromers, which are selectively permeable to bicarbonate and chloride ions. The activation of these channels is triggered by GABA, an endogenous ligand. In neurons, the opening of these channels typically results in chloride ions flowing from the extracellular space into the cytoplasm, leading to a net inhibitory effect [1–3]. GABAARs can be classified into three main categories based on their location and function: (1) postsynaptic receptors, which mediate rapid, point-to-point communication between cells following action potentials; (2) extrasynaptic receptors, characterized by high GABA affinity and a sustained, non-desensitizing ion flow, contributing to tonic inhibition; and (3) perisynaptic receptors, primarily believed to modulate synaptic transmission [4]. Furthermore, the existence of presynaptic GABAARs has also been documented [5–7].
The human and non-human mammalian genomes contain 19 genes encoding for GABAAR subunits, including α1-6, β1-3, γ1-3, δ, ε, π, ρ1-3, and θ subunits. Combined with variations arising from splicing and RNA editing, the potential number of unique GABAAR pentamers is immense, even when considering the assembly rules identified to date [8–16]. The prevailing model suggests that most receptors are composed of two to three β- (or β-like) subunits, one or two α-subunits, and a single ‘odd’ subunit, most commonly γ2, arranged in a counter-clockwise configuration, typically in an α-β-α-γ-β order. However, a significant number of receptor subtypes with unknown or variable native compositions, assembly patterns, and stoichiometries persist [17, 18]. The physiological roles and pharmacological properties of these subtypes exhibit considerable variation, as demonstrated in heterologous expression systems and through in vitro and in vivo studies in rodent models [19–24].
1.2. Benzodiazepines, Z-Drugs, and the Grapefruit Juice Interaction
Benzodiazepines have held a dominant position in the pharmaceutical market for GABAAR-targeting compounds since their introduction by Hoffmann La Roche in the 1960s [25]. They emerged as a safer and more pharmacologically favorable alternative to barbiturates, the previous generation of central nervous system (CNS) depressants acting on GABAARs. Benzodiazepines are heterocyclic molecules defined chemically by an aromatic benzyl ring fused to an unsaturated diazepine ring (Figure 1). Clinically, compounds incorporating a 1,4-diazepine structure are most commonly used. Their effects at low doses are mediated by high-affinity binding sites on GABAARs, located at extracellular interfaces between an α1-3,5 “principal” subunit and a γ1-3 “complementary” subunit [27–33]. Preclinical research and the limited abundance of the γ3 subunit suggest that the four binding sites formed by α1-3,5 in conjunction with γ2 account for the majority of drug effects stemming from these high-affinity interactions. Receptors lacking these high-affinity sites can still be modulated by benzodiazepines at higher (micromolar) concentrations, but they do not display the same low-concentration benzodiazepine effects [34–36]. It’s important to note that this distinction is largely based on data from heterologous expression systems, which may not fully reflect the influence of endogenous GABAAR modulators and their allosteric interactions with benzodiazepine effects [37, 38].
FIGURE 1
Figure 1. Chemical structures of benzodiazepines and Z-drugs, adapted from “New benzodiazepines in Europe review 2021” [26]. (A) Illustrates benzodiazepine scaffolds. (B) Illustrates Z-drug scaffolds, noting that zopiclone is comprised of two enantiomers, not distinguished in this 2D representation.
Benzodiazepines exhibit a wide range of therapeutic actions, including anxiolysis, hypnosis, sedation, muscle relaxation, and anticonvulsant effects [39–42]. Dose-dependent euphoric and amnestic effects have also been reported, potentially contributing to their misuse and illicit use [43]. In the Anatomical Therapeutic Chemical Classification (ATC) System, benzodiazepines are categorized under N03 (antiepileptics) and N05 (psycholeptics), further subdivided into N05B (anxiolytics) and N05C (hypnotics and sedatives). While valuable in treating various conditions like anxiety, epilepsy, and sleep disorders, their broad pharmacological profile often leads to unwanted side effects alongside the desired therapeutic effects. Common side effects include cognitive impairment [44], increased risk of falls and injuries in the elderly [45], sleep architecture disruption [46], sedation, and muscle relaxation [42, 47, 48]. Abrupt cessation after prolonged use can trigger withdrawal symptoms, such as depressive mood, irritability, sleep disturbances, muscle tension, tremors, and even seizures. Benzodiazepine use disorder, while rare, is a serious potential adverse drug reaction. Paradoxical reactions, including talkativeness, restlessness, hyperactivity, agitation, aggressive behavior, and seizures, can occur in less than 1% of patients [46, 49, 50], with juveniles and the elderly being particularly susceptible to adverse effects [39].
Z-drugs—zaleplon, zolpidem, zopiclone, and eszopiclone (see Figure 1)—represent another class of molecules targeting the high-affinity benzodiazepine binding sites on GABAARs. Introduced in the 1990s, they were initially marketed as having lower abuse potential and fewer side effects than benzodiazepines. However, post-launch, adverse event reports associated with Z-drugs have increased. They have been shown to produce euphoric effects comparable to benzodiazepines known for abuse, like lorazepam, alprazolam, and flunitrazepam [51]. Paradoxical reactions similar to those of benzodiazepines have also been reported with Z-drugs [52]. Overall, the side effect profiles of Z-drugs are converging towards those observed with benzodiazepines, although some cognitive measures may show better performance in older populations [53]. For simplicity, we will refer to both benzodiazepines and Z-drugs as “BZ-site ligands”.
The pharmacokinetic and pharmacodynamic profiles of benzodiazepines can be significantly altered by various prescription drugs that interfere with their liver metabolism, particularly through the cytochrome P450 (CYP) system, notably isoenzymes 3A4 and 2C19 [54–57]. This can lead to drug accumulation and severe side effects or, conversely, to therapeutic failure due to accelerated clearance. For example, combining diazepam with drugs like rifampicin or carbamazepine can markedly accelerate its clearance [58–60]. Conversely, hormonal oral contraceptives can reduce clearance and prolong the half-life of several benzodiazepines [61–65]. Notably, grapefruit juice, a common dietary component, can severely inhibit CYP3A4 [66, 67], leading to clinically significant increases in diazepam effects and accumulation. This grapefruit juice benzo interaction is a crucial consideration for patients taking these medications. Furthermore, certain benzodiazepines can interfere with other CYP isoenzyme activities or glucuronidation [68–71].
1.3. Dose Escalation, Non-Medical Use, and Adverse Effect Severity: Controversial Aspects
Between 1996 and 2014, the number of adults in the U.S. filling prescriptions for BZ-site ligands significantly increased, from 8.1 million (4.1%) to 13.5 million (5.6%) [72]. Concurrently, the total quantity dispensed tripled, and overdose deaths involving BZ-site ligands quadrupled from 0.58 to 3.07 per 100,000 adults. Prescription rates for BZ-site ligands show a notable sex disparity, with women receiving prescriptions approximately twice as often as men [73, 74]. Despite this and their widespread clinical use, systematic studies specifically examining sex differences in BZ-site ligand effects are scarce. Existing studies present conflicting evidence regarding sex as a risk factor for substance misuse, with some indicating male sex as a risk factor [75–77] and others suggesting the opposite [78–80]. However, variations in study designs make direct comparisons challenging.
Due to their broad in vivo effects, BZ-site ligands are prominent among commonly misused drugs, categorized as “downers” due to their primarily CNS depressant effects [81]. In the U.S., all benzodiazepines are classified under schedule IV of the “Controlled Substances Act,” indicating a perceived relatively low addictive potential while serving medical needs. Nonetheless, both benzodiazepines and Z-drugs can induce physical and psychological dependence even after short periods, necessitating the principle of “as short as possible, as long as needed” for treatment regimens. Drug tolerance, requiring higher doses for the same effect, is another concern associated with BZ-site ligands.
While generally considered a safe class of compounds, benzodiazepines and Z-drugs can pose problems in long-term treatment and illicit drug use. Benzodiazepine overdose alone carries a comparatively lower risk of fatal outcomes than other depressants like barbiturates, but the risk is not negligible [82–84]. Emergency department visits involving BZ-site ligands nearly doubled between 2005 and 2011, according to the DAWN (Drug Abuse Warning Network) report [84]. Benzodiazepine users face a higher risk of serious outcomes during emergency visits, further amplified by combining benzodiazepines with alcohol or opioids [84]. Benzodiazepines have also been linked to an approximate doubling of motor vehicle accident risk [85], with similar effects reported for zopiclone [86, 87].
BZ-site ligands are often used in combination with other drugs, rather than as primary drugs of abuse [43]. Benzodiazepines with rapid onset, such as diazepam and alprazolam, can produce euphoric effects at higher concentrations and are sometimes combined with methadone to enhance its mood-altering effects [43]. Cocaine and other stimulant users may use benzodiazepines to mitigate side effects or for “coming down” [43]. Analysis of over 1,200 oxycodone-related drug abuse deaths revealed a high prevalence of diazepam co-abuse [84], a finding corroborated by other sources [88]. The combination of alcohol and BZ-site ligands is particularly concerning given alcohol’s social acceptance and easy accessibility [89]. Both alcohol and BZ-site ligands are CNS depressants, leading to compounded effects when taken together. Individuals with alcohol use disorder (AUD) may experience a stronger psychoactive effect from benzodiazepines. Those with a family history of AUD may also exhibit altered sensitivity and effects to alprazolam [90–93]. While studies document drug-alcohol interactions and adverse outcomes linked to BZ-site ligands, systematic comparisons between individual drugs are lacking. Consequently, it remains unclear whether certain substances are more problematic in different forms of medical and non-medical use, despite anecdotal evidence suggesting specific compounds are particularly prone to misuse, such as in date rape scenarios or triggering negative experiences.
Despite these concerns, Bz-site ligands remain indispensable in daily clinical practice [94]. Attempts to restrict their use through stricter prescription regulations have paradoxically led to increases in overdose emergencies involving drugs with less favorable safety profiles [95–97]. Furthermore, acute sedation or respiratory depression from excessive Bz-site ligand doses can be readily reversed with intravenous flumazenil in clinical settings. Therefore, enhancing awareness and understanding of the mechanisms behind unwanted and potentially dangerous benzodiazepine effects is crucial for rational clinical decision-making and for developing safer drugs with similar therapeutic benefits.
1.4. FAERS and the Utility of Pharmacovigilance Data
Comprehensive, controlled studies and head-to-head comparisons of individual benzodiazepines and Z-drugs according to current scientific standards are lacking for most indications, partly because many of these substances are no longer patent-protected.
Pharmacovigilance, a rapidly evolving scientific field, is dedicated to detecting, assessing, understanding, and preventing drug-related problems and adverse events, encompassing all activities aimed at improving drug safety [98]. The U.S. Food and Drug Administration (FDA) maintains the FDA Adverse Event Reporting System (FAERS), a federal database for collecting real-world, post-marketing observations of drug-related adverse events. FAERS includes reports of adverse events, medication errors, and product quality complaints submitted by healthcare professionals, patients, and the public [99]. Pharmacovigilance analysis typically uses metrics derived from this data to assess the association between a drug and an adverse event. Higher values for these metrics indicate a stronger association. Post-marketing pharmacovigilance data like FAERS can provide valuable signals for adverse reactions not detected during the initial drug approval process. Due to inherent limitations of real-world data, such as incomplete dosage information, co-medication details, and demographic data, specialized data analysis methods have been developed [100–102]. While a strong signal suggests an association between a drug and an outcome, it does not establish causation. Pharmacovigilance data is not suited for comparative pharmacology [101]. However, it serves as the only available data source for generating hypotheses based on large-scale, real-world observations across a wide range of drugs.
FAERS uses “preferred terms” to report adverse events, categorized using the MedDRA (Medical Dictionary for Regulatory Activities) terminology (Figure 2). MedDRA is a hierarchical medical terminology system that allows analysis of FAERS reports at different levels. The five levels are System Organ Class (SOC), High Level Group Term (HLGT), High Level Term (HLT), Preferred Term (PT), and Lowest Level Term (LLT) [103]. Figure 2A illustrates the MedDRA hierarchy and the levels used in this study.
FIGURE 2
Figure 2. MedDRA dictionary system overview and IC025 usage. (A) The five MedDRA hierarchy levels with examples; red: System Organ Class (SOC); blue: High Level Group (HLG), Cyan: Preferred Term (PT). Arrows indicate direction from higher to lower levels, dashed arrows show analysis surpassing High Level Terms. (B) IC025 values and calculations. IC025 reflects drug-AE association; cIC025 is the sum of IC025s for a drug within an HLG; ScIC025 is the sum of cumulative IC025s for an HLG or SOC.
Given the widespread prescription of Bz-site ligands (benzodiazepines and Z-drugs), a large volume of reports are associated with their use. However, a comprehensive comparison of individual compounds within this class based on these reports has been lacking. This study employs disproportionality analysis, providing mathematically defined parameters for the strength of association signals [104]. Specifically, the information component (IC) value, and in particular IC025 (the lower end of the 95% credibility interval [100]), quantify the strength of the dependency between a specific drug and a reported adverse event (Figure 2). Analyzing a large FAERS dataset, we aimed to generate individual drug profiles and identify tendencies of drug heterogeneity, some of which are corroborated by clinical study data.
Molecular mechanisms underlying drug heterogeneity may arise from both off-target and on-target effects. While exhaustive investigation of every potential off-target effect is impractical, on-target heterogeneity can be explored through structural data, providing hypotheses for mechanisms driving the diverse effects of these drugs. These mechanisms largely stem from multiple binding sites and their cooperative interactions at various receptor subtypes. Therefore, compounds with informative FAERS records were also examined for chemical features influencing pharmacodynamics, aiming to identify drug properties associated with specific unwanted effects. A brief overview of binding site heterogeneity within the GABAAR family is also provided.
2. Results
We mined the publicly accessible FAERS (FDA Adverse Event Reporting System) dataset from Khaleel et al. [105] to establish pharmacovigilance profiles for each drug and sex, focusing on Bz-site ligands with sufficient data. Datasets for female and male reports were analyzed separately, similar to the approach used by Drug Central [106]. The analysis workflow is depicted in Figure 3 (and detailed in the Methods section). The complete dataset comprised 170,565,117 drug-adverse event combinations, with 100,085,277 for females and 59,680,210 for males. These reports were filtered for 173 benzodiazepines and Z-drugs from our predefined list (see Methods and Supplementary Item 1). This filtering resulted in 2,701,733 female and 1,447,028 male drug-adverse event combinations, representing 44 benzodiazepines with FAERS entries, designated as data pool 1 (Figure 3; Supplementary Table S1). Disproportionality analysis was then performed to identify drug-adverse event associations. Data pool 2 inclusion criteria were PRR > 2 and IC025 > 0 [107–110], along with a minimum of five reports per drug-AE pair. Supplementary Item 2 details the composition of pools 1 and 2 (Figure 3) regarding total reports per drug. The IC025 value served as the primary criterion for subsequent data filtering, as recommended by the UMC (Uppsala Monitoring Center) [111], with higher IC025 values indicating stronger signals. Of the initial 44 drugs, 39 met the criteria for further analysis (pool 2). Raw data is available in Supplementary Items 3, 4 (Excel format).
FIGURE 3
Figure 3. FAERS dataset analysis pipeline. Green boxes: results after filtering steps (blue boxes). Data for MedDRA category analysis is “pool 1” and “pool 2”.
Pool 2 records (post-disproportionality analysis), containing only drug-AE associations, were analyzed using MedDRA categories. Pool 1 data was occasionally used for comparison. The nomenclature used for pool 2 data is as follows (Figure 2B): IC025 represents the value for a specific drug and individual adverse effect combination, with only positive values from pool 1 used for aggregate value calculations. Cumulative IC025 (cIC025) is the sum of all positive IC025 values for a specific drug within a category (HLG or SOC). Summative cumulative IC025 (ScIC025) is the sum of cIC025 values for all drugs within a group (HLG or SOC) (Figure 2B).
2.1. Overview Across all SOCs
Across all System Organ Classes (SOCs), associations were identified for all 39 drugs in 27 SOCs (Figure 3). To gain an overview, summed cumulative IC025 (ScIC025) values per SOC were calculated and are displayed in Figure 4A. As expected, the highest summed cumulative signals were observed for “nervous system disorders” and “psychiatric disorders,” combined as the “neuropsychiatric” group. The high ScIC025 for “injury, poisoning and procedural complications” (Figure 4A) is likely due to the widespread non-medical use of Bz-site ligands, followed closely by “investigations”. The top four SOCs were further broken down into their contributing HLGs (Supplementary Figures S1–S3).
FIGURE 4
Figure 4. AE association distribution across organ system classes: (A) Summed cumulative IC025 (ScIC025) values per SOC for drugs with positive IC025, separated by sex (female: light red, male: light blue). (B) Pie charts for top 10 drugs per sex (total ≥10) with highest cIC025 contribution to top four SOCs by summed cIC025 (ScI025). Segment size corresponds to cIC025 contribution (percentage). Outer circles: males, inner: females, drugs sorted by female ScI025 rank (clockwise, descending).
For each of the top four SOCs, drugs were ranked by ScIC025 to assess gross contributions. The top 10+ drugs for each SOC are shown in Figure 4B, with >10 drugs shown due to sex-based differences in top 10 lists. Each SOC exhibits a unique drug ranking, with different top-ranked drugs. Clobazam ranks highest in “nervous system disorders” for both sexes and is in the top 10 for the other three SOCs. Clonazepam is top-ranked in “psychiatric disorders” and is also in the top 10 across all four SOCs. The relative contribution of each drug to each SOC varies between sexes. Midazolam, used as a procedural anesthetic, is top-ranked in “injury, poisoning and procedural complications” but not in the top 10 for “psychiatric disorders”. “Investigations” is a heterogeneous SOC, reflecting parameter changes across systems, not a single organ system. Striking cIC025 differences between sexes are seen for temazepam and nitrazepam in this SOC. Closer examination (Supplementary Figure S3) reveals high female ScIC025 in cardiac investigations, consistent with higher female ScIC025 in “cardiac disorders” (Figure 4A).
Results from the top four SOCs suggest heterogeneous side effect patterns across compounds. Some drugs appear only in the top 10 of specific SOCs (e.g., brotizolam in “investigations”), while others dominate multiple or all SOCs. To further investigate drug heterogeneity, we focused on the top two SOCs, merging them into “neuropsychiatric reports” (Figure 5). Supplementary Figures S1, S2 provide separate details for “nervous system disorders” and “psychiatric disorders”.
FIGURE 5
Figure 5. AE association distribution across “nervous system disorders” and “psychiatric disorders”: Polar bar charts show Higher Level Groups (HLGs) with summed cumulative IC025 ≥300. cIC025 contribution per compound is shown as patterned segments (legend). Drugs are sorted by total ScIC025 contribution within the SOC per sex (high to low, outer rims). Table below: nine largest contributions from top two SOCs. Green: nervous system disorders; blue: psychiatric disorders.
2.2. Analysis of Neuropsychiatric AEs
Further analysis focused on neuropsychiatric SOCs. The largest contributing higher level groups (HLGs) within these SOCs are shown in Figure 5. These include “neurological/psychiatric disorders not elsewhere classifiable (nec)”, two groups related to movement/motor disturbances, sleep-related signs and symptoms, anxiety signs, seizures, suicidal and self-injurious behaviors, and disturbances in thinking and perception. ScIC025 per HLG shows minor sex differences (Figure 5). Next, drug contributions (cIC025) to each HLG were examined.
Consistent with their large ScIC025 contributions to overall nervous system and psychiatric SOCs, clobazam, clonazepam, and diazepam also show considerable cIC025 values in individual neuropsychiatric HLGs. However, heterogeneity is also evident at this level. Clobazam shows a strong association with the “seizures” HLG, consistent with previous findings [112], but a smaller signal for “anxiety disorders and symptoms”. In movement and muscle symptom groups (Figure 5B), drugs exhibit varying association strengths, with clobazam showing a stronger signal in males.
2.3. Neurological Disorders Not Elsewhere Classified
The HLG “neurological disorders nec” (Figure 5A) contributes almost half of the ScIC025 from the “nervous system disorders” SOC, with summed cumulative IC025 values of approximately 900 in females and 800 in males. This HLG was analyzed in detail at the individual adverse event association level (IC025 values) (Figure 6).
FIGURE 6
Figure 6. AE association distribution across “neurological disorders nec” HLG: Polar bar charts show adverse events with signal in this HLG as cumulative IC025 across all contributing drugs (cutoff: cIC025 > 30). Full dataset: Supplementary Items 3, 4. Drugs identified in list below graph. Panel B enlarges panel A (cIC025 scale 0–20).
The main contributing AEs to this HLG are sedation and over-sedation signs, including sedation, somnolence, sopor, and coma, indicating varying degrees of CNS depression (Figure 6). The second largest group, encompassing agitation, restlessness, and logorrhea, reflects paradoxical responses (Figure 6). Interestingly, the cumulative signal for logorrhea differs significantly between sexes, with drug specificity. Logorrhea is prominent in the male dataset (IC025 = 5.3) but absent for females, due to low report numbers (Supplementary Items 3, 4). A large portion of cIC025 in neuropsychiatric groups reflects expected sedation/over-sedation signs and paradoxical reactions.
2.4. Highest Ranked Psychiatric HLGs
To explore contributions beyond sedation and paradoxical responses, we analyzed the top four psychiatric HLGs by summed cumulative IC025 (Figure 3A) individually (Figure 7). These include “psychiatric disorders nec”, “sleep disorders and disturbances”, “anxiety disorders and symptoms”, and “suicidal and self-injurious behavior”. Sleep-related disturbances, present in both psychiatric and nervous system disorder SOCs, were merged prior to analysis (Figure 2; Supplementary Figure S4). Individual AEs contributing to each group are in Supplementary Figures S4–S7 and Supplementary Table S2.
FIGURE 7
Figure 7. Detailed analysis of top four psychiatric HLGs. (A) Scatter plot (left axis): cumulative IC025 per drug in each HLG, highest cIC025 drugs labeled. Bar graphs (right axis): report percentage of all drugs within HLG vs. all drugs in all HLGs (blue: male, red: female). (B) Percentage of reports per drug for each HLG relative to total reports per drug (color legend). Raw data: Supplementary Items 2–5 and Supplementary Table S2.
The HLG “psychiatric disorders nec” has a higher fraction of total reports for males vs. females, while the other three HLGs show relatively balanced values (Figure 7A). Diazepam has the highest cumulative signals for both sexes in “psychiatric disorders nec” and “suicidal and self-injurious behavior”. Diazepam, alprazolam, and clonazepam are among the top five drugs for both sexes (Supplementary Table S2). Figure 7B shows the fraction of reports within each HLG from panel A for each drug (normalized to 100% per drug). Nordazepam contributes significantly to the male “psychiatric disorders nec” dataset, with >13% of its associated AEs from these four HLGs (Figure 7B). “Psychiatric disorders not elsewhere classifiable” mainly includes AE associations related to abuse and withdrawal (Supplementary Figure S5). Nordazepam generally shows a strong signal for “drug abuse,” “substance abuse,” and related AEs. Oxazepam and lormetazepam also have high IC025 values for most abuse/addiction-related AEs compared to loprazolam and flunitrazepam, which have weak or no signals. Alprazolam is uniquely associated with all investigated AEs in this category (Supplementary Figure S8).
Clonazepam has a strong cIC025 in “psychiatric disorders nec” and ranks in the top 5 for “sleep disorders and disturbances,” “suicidal and self-injurious behaviors nec,” and “anxiety disorders and symptoms”. These four AE groups account for nearly 10% of all AEs associated with clonazepam in males (Figure 7B).
In the merged HLG “sleep disorders and disturbances”, eszopiclone and zolpidem show the strongest signals. Eszopiclone has an exceptionally high fraction of reports in this HLG (Figure 7B). At the individual AE level, eszopiclone is primarily associated with insomnia signs, while zolpidem is linked to sleep disturbances like somnambulism and sleep-related eating issues (Supplementary Item 5). In sleep disorders and disturbances, Z-drugs (zolpidem, eszopiclone, zaleplon) are among the top five ranked drugs (Supplementary Table S2).
For the HLG “anxiety disorders and symptoms”, clonazepam is top-ranked for both sexes. Normalized report numbers show sex differences for some drugs (Supplementary Figure S7. Eszopiclone has a high cumulative signal for females but not males, while midazolam, tetrazepam, and zaleplon show stronger associations for males. The top AE in this HLG, agitation (Supplementary Figure S7), is also in the “neurological disorders nec” group and not anxiety-specific. However, agoraphobia and panic signs contribute significantly to clonazepam’s signal strength in both sexes.
These examples highlight varying drug association strengths with different symptom groups and sex-specific differences. While requiring careful interpretation, the data strongly suggests significant compound heterogeneity, indicating unexpected differences in side effect profiles among BZ-site ligands.
2.5. Sex Differences
Given the observed sex-specific signals, we analyzed neuropsychiatric SOC data for sex differences per drug. Neuropsychiatric data was visualized in a scatter plot of sex-specific IC025 values (Figure 8A). Drug/AE pairs with a ≥2:1 IC025 ratio were highlighted to identify sex-biased associations. Neuropsychiatric drug/AE events exclusively in one sex are shown as points on the axes (Figure 8A). AEs with the highest cumulative IC025 from drugs with a >2:1 signal ratio were analyzed for each sex (top 20 each, Figures 8B,C). Drug-AE pairs with signals only in one sex were also analyzed (top 20, Supplementary Figures S9, S10). The top 10 drugs contributing most to adverse events in the same dataset (ratio > 2:1) were also identified (Figure 8D).
FIGURE 8
Figure 8. Sex differences in neuropsychiatric AE signal strength: (A) Scatter plot of AE/drug pairs, strong colors for ≥2:1 IC025 ratio, pale for <2:1. (B,C) Top 20 neuropsychiatric AEs (male: B, female: C) from data in panel A (≥2:1 ratio), drug contributions color-coded. (D) Top 10 contributing drugs from male and female data (≥2:1 ratio).
Several psychiatric AEs show sex-biased signals. Paradoxical response signs like logorrhea, mania, restlessness, and hypomania are more strongly associated with males but absent from top-ranked female AEs (>2:1 dataset). Nordazepam has a notably strong male signal for paradoxical responses and various neuropsychiatric signs (Supplementary Figure S11). Logorrhea is exclusively associated with male reports for nordazepam, tetrazepam, and oxazepam (Supplementary Figure S9). For female:male ratios ≥2:1, top-ranked AEs are mixed, including unspecific neuropsychiatric changes. Anterograde amnesia is strongly associated with clorazepate and lormetazepam in females, with no male association (Supplementary Figure S10). Intrigued by this, we investigated amnestic effects and confirmed strong associations for these two compounds, absent in males (Supplementary Figure S12). Zolpidem shows amnestic AE associations in both sexes, more so than classical benzodiazepines.
Zolpidem is among the few compounds with stronger male:female signals (>2:1), while bromazepam and lormetazepam show the opposite. Most compounds with stronger male associations have more female reports in pool 2, suggesting a specific phenomenon, not just reporting bias. These sex-specific AE profiles from FAERS data warrant further pharmacoepidemiological and clinical follow-up.
2.6. Physico-Chemical Descriptors
Pharmacovigilance data strongly suggests that Bz-site ligand profiles differ significantly in in vivo human effects, and safety profiles for some drugs might be incomplete. Molecular drivers of this heterogeneity are complex and deserve consideration to align pharmacovigilance findings with testable hypotheses for future research.
To understand molecular similarity and heterogeneity patterns at the drug level, we evaluated compounds based on physico-chemical and 3D properties. Integrating pharmacophore models and fingerprints with pharmacovigilance data can reveal previously unknown safety concerns related to drug scaffolds, enabling measures to improve drug safety, such as avoiding specific moieties in drug development or restricting certain drugs in specific patient groups. Systematic analysis of drug similarities and differences can also identify unexpected structure-activity landscape cliffs, informing bed-to-bench considerations for drug development.
Establishing structure-activity relationships for GABAAR-targeting molecules has been challenging due to small chemical scaffold changes causing unexpected heterogeneity in structure-activity landscapes from heterologous expression systems and in vivo outcomes [113, 114]. Here, we approached heterogeneity from a pharmacovigilance perspective and related findings to ligand-based approaches. We used ligand fingerprints (physico-chemical parameters) and pharmacophore models (3D shape/functional groups) to describe molecule properties and cluster compounds (Figure 9).
FIGURE 9
Figure 9. Clustering of 39 benzodiazepines and Z-drugs (pool 2, Figure 2) by physico-chemical parameters and 3D properties. Upper panel: ligand fingerprint analysis (physico-chemical descriptors). Lower panel: pharmacophore approach (3D shape/functional groups). Stereoisomers denoted by “_1” or “_2”. Values from vROCS® (Supplementary Item 6). Connecting lines: compounds grouped by both methods. Dashed lines: selected clustering differences.
Ligand fingerprint analysis revealed a more complex molecular landscape than expected from 2D/3D-structure similarity, consistent with pharmacophore analysis. Pharmacophore analysis showed distinct clusters for “xazolam” compounds (cloxazolam, ketazolam, mexazolam, oxazolam), triazolo-compounds, and 1,4-benzodiazepines (diazepam metabolites) (Figure 9, lower panel). Ligand fingerprint analysis generally agreed but showed unexpected clusters, like zaleplon with lormetazepam (Figure 9, upper panel), structurally dissimilar in 3D. Zolpidem clustered with midazolam and clotiazepam despite different 3D shapes. Ligand fingerprint analysis, based on physico-chemical properties, provided a more nuanced, less intuitive picture complementing 3D pharmacophore results.
2.7. Complexity of On-Target Effects
FAERS signals suggesting compound heterogeneity are not surprising given past observations of compound promiscuity at single GABAARs and the existence of multiple homologous GABAAR subtypes [17, 115, 116]. Recent structural data confirms non-canonical binding sites and differential binding modes, providing structural correlates for compound-specific effects. The current structural evidence is summarized in Figure 10.
FIGURE 10
Figure 10. Benzodiazepine and zolpidem binding sites. (A) Extracellular high affinity binding site (BS) shared by BZDs and Z-drugs, requiring α1,2,3, or 5 subunit (+) and γ1-3 (−) interface. Light gray: α1+; dark gray: γ2-; red surface: ligands bound in structures (table). (B) Diazepam and zolpidem low affinity binding sites in transmembrane domain (TMD). Diazepam in α2+/β2- and β2+/α1- pockets. Zolpidem in β2+/α1- (8DD2 [117]). Light gray: β2+/γ2+; dark gray: α1−/β2-; ligand surfaces: blue/cyan. Insert: ligand superposition, color-coded binding modes.
Structural evidence highlights: (i) Common core compounds can have different high-affinity site binding modes (flumazenil vs. diazepam/alprazolam [118–120]). (ii) Non-canonical sites predicted by mutagenesis are confirmed and expanded by structural evidence [121, 122], showing distinct binding modes (Figure 9B) [120]. (iii) Biochemical evidence for non-canonical sites (ECD β2+/γ2-interfaces [123]) is supported by receptors lacking alpha subunits [17].
2.8. Chemical Similarity and Pharmacological Trends: A Complex Relationship
Drug-protein interactions for chemically similar ligands form structure-activity landscapes with “smooth” and “rugged” features. Chemically clustered compounds are more likely to have overlapping pharmacological profiles to some extent. Intrigued by dissimilar AE patterns for zopiclone and eszopiclone, we analyzed them with chemically similar brotizolam and etizolam (Figure 11). Eszopiclone shows exceptionally high IC025 and fractional report share in “sleep disorders and disturbances”, while zopiclone lacks strong sleep disturbance associations. This is intriguing as zopiclone is racemic eszopiclone and less affine/active R-enantiomer [124, 125].
FIGURE 11
Figure 11. Top 10 neuropsychiatric AEs of four drugs examining chemical similarity patterns. (A) Eszopiclone and zopiclone: Top 10 AEs for eszopiclone (IC025) and zopiclone (values added) (left); top 10 AEs for zopiclone (IC025) and eszopiclone (values added) (right). Y-axis: IC025 (positive for association). Negative values: pool 1 data. (B) Brotizolam and etizolam: Top 10 AEs for brotizolam (IC025) and etizolam (values added) (left); top 10 AEs for etizolam (IC025) and brotizolam (values added) (right). Y-axis: IC025 (positive for association). Negative values: pool 1 data.
Eszopiclone and zopiclone top 10 AEs show no overlap. Eszopiclone’s strong AE associations are sleep and sleep-related disturbances, affecting both sexes similarly. Zopiclone shows mixed positive/negative associations for sleep AEs. Zopiclone’s top 10 AEs cover a broad range of neuropsychiatric phenomena, some absent for eszopiclone. Dystonic tremor is male-specific, muscle spasticity shows negative male and positive female associations. Intentional self-injury is strongly associated with zopiclone; eszopiclone’s negative IC025 suggests this is zopiclone-specific. This suggests R-zopiclone is not just lower affinity [124] but has specific, dominant side-effects and may overcome eszopiclone’s paradoxical responses on sleep side effects. Etizolam and brotizolam, closely clustered chemically and by pharmacophore, show highly similar FAERS profiles, as expected.
3. Discussion
FAERS data limitations require pharmacoepidemiological studies to validate associations. We examined selected drug-AE combinations with strong signals in our analysis against literature and databases like SIDER [126] (product information derived from clinical trials). Results are in Table 1.
TABLE 1
Table 1. FAERS signals and other evidence for selected drug-AE event pairs: SIDER database: “frequent”/“common” (+), others (+), not mentioned (-), n/a: drug absent.
This evidence compilation, while not exhaustive, shows strong FAERS associations are often confirmed in systematic studies.
This study challenges the notion of benzodiazepines and Z-drugs as interchangeable, except for pharmacokinetic differences, contrasting with anecdotal evidence and early preclinical literature [133]. Comprehensive FAERS analysis of BZ-site ligand AEs (2004Q1-2021Q3) reveals diverse AE portfolios per compound, varying significantly between drugs. This is partly reflected in product information and literature, but systematic data and safety profiles are limited. Pharmacovigilance data is a valuable source of human observations for identifying individual compound profiles.
We focused on neuropsychiatric AEs, as BZ-site ligands are mostly psychotropic, except for some antiepileptics. In neuropsychiatric MedDRA terms, only 11 of 39 drugs account for >58% of total neuropsychiatric ScIC025 (Supplementary Figure S13), suggesting differing neuropsychiatric AE propensities. This isn’t primarily due to prescription bias, as drugs with low (nordazepam) and high (diazepam) report numbers, and highly prescribed drugs with low neuropsychiatric signals (triazolam, eszopiclone) were all found (Figure 5, Supplementary Figures S1, S2).
Outside neuropsychiatric SOCs, “investigations” HLG showed high cIC025. Female records showed particularly strong changes in electrocardiogram parameters and blood pressure compared to males (Supplementary Figures S3, S14, S15).
Over-sedation and related effects, and paradoxical responses (“neurological disorders nec” group, Figure 6) dominate total neuropsychiatric ScIC025. Paradoxical responses may be drug-specific and should be quantified per drug, not class-wide. Pharmacovigilance data shows drug specificity, consistent with preclinical work. The striking eszopiclone (hypnotic) association with insomnia and sleep disturbance signs, nearly absent in zopiclone data (Figures 7, 11), is notable.
In psychiatric HLGs, beyond over-sedation effects, anxiety symptoms and self-injurious behaviors show strong signals. Another zopiclone/eszopiclone difference: self-injurious behaviors are strongly associated only with zopiclone (Figure 11). FAERS data suggests substantial drug heterogeneity, affecting both medically used BZ-site ligands and non-medical use issues. We extracted IC025 values per drug for dependence, drug abuse, and withdrawal symptoms (Supplementary Figure S8). While requiring careful interpretation, this data reveals drug heterogeneity warranting further investigation.
Non-medical drug use can be recreational (downers, enhancing other drugs) or illicit (date rape). Amnestic effects are relevant, with anterograde amnesia associated with female records exclusively for lormetazepam and clorazepate (IC025 > 3, Supplementary Figures S10, S12). Only 17 of 39 drugs were associated with any amnestic effects; zaleplon, for example, showed none. Drugs show heterogeneity regarding properties compatible with date rape drug abuse. While FAERS data has limitations, these hints of drug heterogeneity should prompt further systematic investigations.
This study highlights pronounced drug differences in AE profiles for many compounds. Given women receive twice as many benzodiazepine prescriptions in the US, a similar AE report ratio in FAERS might be expected absent sex-specific factors. However, data interpretation is limited by: (1) Sex-specific reporting bias; (2) Sex-biased prescription patterns for certain benzodiazepines (e.g., anxiety disorders in women); (3) Higher male non-prescription benzodiazepine abuse; (4) Missing drug and sex-specific prescription numbers in FAERS; (5) Reports from different indications and dosages; (6) Illicit use reports lacking prescription and often involving polydrug use; (7) Other confounding factors.
However, preclinical and in vitro research suggests mechanisms for sex differences beyond data bias: Supra-additive effects with endogenous cannabinoids [37] and (neuro-) steroids [134] have been observed for benzodiazepines in vitro. Native GABAARs may contain endogenous allopregnanolone before benzodiazepine addition [38].
We attempted to correlate molecular properties with FAERS-derived drug profiles but found little evidence for correlations between chemical properties and pharmacovigilance fingerprints. The vastly different eszopiclone/zopiclone FAERS profiles demonstrate ligand-based approach limitations. Our data suggests both: Highly similar compounds can share drug properties (etizolam/brotizolam, Figure 11), and steep structure-activity landscape cliffs can occur (S- and R-zopiclone). The latter is structurally correlated with multiple binding sites with varying affinity and efficacy. The broad range of in vivo effects observed for benzodiazepine site ligands is likely mediated by multiple distinct receptor subtypes with unique binding sites [135]. Rodent subtype-specific pharmacology research has had limited translational success [113], partly due to neuronal cell type transcriptome differences and regio-specific subunit expression discrepancies between animals and humans [136]. Evidence suggests that both “high affinity” and additional binding sites (shared with general anesthetics) contribute to in vivo effects. Differences in compound ability to utilize these sites lead to specific pharmacodynamic profiles. Individual BZ-site ligands vary in their tendency to occupy additional sites beyond canonical high-affinity sites [120, 122], likely impacting in vivo effects due to near-complete isoform-indifference in low-affinity sites [137]. Recent structural findings offer an updated view of known and putative allosteric sites mediating pharmacological effects (Figure 10).
Furthermore, compound binding and effect heterogeneity at canonical sites is understudied. γ1 and γ3 pharmacology is incomplete, and high γ1 subunit expression in the human limbic system may explain specific drug effects for substances acting on γ1-containing receptors [133]. Modulatory efficacy, ranging from strong GABA enhancing (PAM) to diminishing (NAM) effects [133, 138–140], adds complexity. Data on modulatory effects in major receptor subtypes is lacking for most approved benzodiazepines and Z-drugs, with PAM effects assumed (except for flumazenil). Beyond diverse allosteric sites on GABAARs, off-target effects may also drive individual drug effects, although CNS-target assays suggest fewer off-targets for these compounds compared to other CNS therapeutics.
In summary, this study provides insights into BZD and Z-drug pharmacological properties, informing clinical decision-making and drug development. FAERS dataset analysis and ligand fingerprint/pharmacophore analyses reveal BZ-site ligand heterogeneity, aiding in identifying therapeutic uses and adverse effects and shaping future clinical studies. FAERS profiles suggest sex-specific side effects.
Our analysis yielded strong drug-AE associations. While pharmacovigilance data doesn’t confirm causality, our findings may stimulate follow-up research and potentially adapt prescription practices for sex-specific care. If clinical studies confirm FAERS associations, product information and legal classifications could be adjusted to reflect increased risks of unwanted effects, including specific warnings.
4. Materials and methods
4.1. Data Mining
Four public sources were used to compile a list of benzodiazepines and Z-drugs: Drugbank [141], Wikipedia [142, 143], and Wikidata [144]. Different data extraction methods were used for each source. For Drugbank, a Python script filtered compounds linked to each GABAAR subunit, manually screened for benzodiazepines and Z-drugs. SPARQL was used to search Wikidata for drugs related to any of the 19 subunits. Benzodiazepines were also obtained from two Wikipedia pages [142, 143]. The final list included 173 benzodiazepines and Z-drugs.
4.2. FAERS Analysis
Pharmacovigilance analysis of benzodiazepines and Z-drugs used a FAERS (FDA Adverse Event Reporting System) dataset from Khaleel et al. [105], covering Q1 2004 to Q3 2021. The dataset was initially split into male and female subsets, removing records with unknown sex. Female reports concerning offspring occurrences were manually curated and removed if needed. Disproportionality analysis was performed for each drug-adverse event pair in both datasets.
Disproportionality analysis assessed the association strength between drug use and adverse events (AEs) [100–102]. For each drug-adverse event pair, information component (IC), 95% confidence interval of IC (IC025) [100, 102], proportional reporting ratio (PRR), and reporting odds ratio (ROR) were calculated [101, 105]. IC assesses the likelihood of true values within a range. The Uppsala Monitoring Centre developed and validated this Bayesian neural network approach for IC [111, 145], representing the base 2 logarithm of observed/expected ratios, used for WHO database analysis [100, 102, 146]. IC025 has been used as a benchmark for positive drug-adverse event connections [102, 104, 147–151]. PRR measure with IC025 and ≥5 observations were used for affirmative signals [107–110], as in this work. Relevant records for 173 benzodiazepines and Z-drugs were extracted, examining System Organ Classes (SOCs) and Higher Level Groups (HLGs) within MedDRA (version 22.1). Only adverse events identifiable in MedDRA version 22.1 were considered. Individual AE levels were analyzed with filters. Cumulative IC025 values and relative/absolute report numbers were obtained as sums from filtered records.
a = Reports of the drug of interest with the adverse event of interest. b = Reports of all other drugs with the adverse event of interest. c = Total drug reports of all other adverse events. d = Total reports of all other drugs with all other adverse events.
Proportional Reporting Ratio(PRR)=a(a+c)b(b+d)
Information component(IC)=log2a+0.5aexp+0.5
aexp=(a+b)∗(a+c)(a+b+c+d)
IC025=IC−3.3∗(a+0.5)−12−2∗(a+0.5)−32
Reporting Odds Ratio(ROR)=acbd
Python libraries pandas 2.8.2 and numpy 1.22.4 with python 3.10 were used for calculations. Figures were generated with plotly 5.13 with python 3.10.
Calculation for Figure 4B:
Calculation:
Drug percentage=cumulativeIC025ofaspecific drug inSOCsumof cumulativeIC025ofalldrugs inSOC∗100.
Calculation for Figure 7B:
Calculation:
Reports percentage=total reports of drug inHLGtotal reports of drug∗100
4.3. Ligand Based Methods and Plot Creation
4.3.1. Fingerprints
Python 3.10 and PyBioMed 1.0 library (moe, ghosecrippenfingerprint, cats2d, connectivity, topology) were used to calculate molecular fingerprints from 39 drug sdf files. Principal component analysis (PCA) was performed using scikit-learn 1.0.2 library. Plotly version 5.13 created dendrograms for graphical representation.
4.3.2. 3D-structure similarity analysis
3D structures of drugs were retrieved from PubChem [152] as SD-files, merged into a single file. Flipper software [153] (version 3.1.1.2) enumerated stereoisomers (settings: enumEZ true-enumNitrogen false-enumRS true-enumSpecifiedStereo true-warts true). Stereoisomers are named with an underscore and number suffix. Omega software [153, 154] generated conformer ensembles (settings: -ewindow 20.0, -rms 0.25). Conformers were stored as a single SD-file, used for shape-based drug similarity calculations using ROCS program [154]. Shape similarity values for all drug stereoisomer pairs were calculated using the multi-conformer SD-file as input for query (-query) and database (-dbase) structures. ROCS parameters used defaults, except: -mcquery false, -scdbase true, -report one. A Python script (‘report_to_dist_matrix.py) converted ROCS report file to ‘shape distance’ matrices based on ColorTanimoto, ShapeTanimoto, and TanimotoCombo scores. Only TanimotoCombo was used further. Shape distance Dij for drug stereoisomer pair ij was calculated as:
Dij = 1 – Sij/Smax
Smax is the maximum similarity score (1.0 for Color- and ShapeTanimoto, 2.0 for TanimotoCombo).
ROCS combo scores were used for further analysis. Dendrograms were created using python 3.10 with plotly 5.13.0.
4.4. Analysis of Structural Data
PDB was mined for GABAA receptor structures with analyzed BZ-site ligands. Structures [8DD2, 6X3X [118], 6HUO [119], 6X3U [120]] were superposed and rendered with Schrödinger/Maestro Version 13.1.141.
Data availability statement
The original contributions are publicly available at: https://github.com/FilipKon/Drug_Diversity_FAERS.
Author contributions
FK and TS performed experiments. FV, FK, ID, MW, MK, TS, and ME contributed to writing. FV, ME, MW, ID, and FK shaped manuscript scope and literature research. ID and MW chose drugs. FV and ME designed experiments. ME provided funding and supervised. All authors contributed and approved the submitted version.
Funding
The authors acknowledge funding from the Austrian Science Fund (MolTag doctoral programme FWF W1232) and the European Community (NeuroDeRisk project, Innovative Medicines Initiative 2 Joint Undertaking, grant agreement No. 821528). This Joint Undertaking receives support from the European Union’s Horizon 2020 and EFPIA.
Acknowledgments
Discussions with Jure Fabjan on data analysis and visualization are acknowledged.
Conflict of interest
The authors declare no commercial or financial conflicts of interest.
Publisher’s note
All claims are solely those of the authors and do not represent affiliated organizations, publisher, editors, or reviewers. Product evaluations or claims are not guaranteed or endorsed by the publisher.
Supplementary material
Supplementary material is available online at: https://www.frontiersin.org/articles/10.3389/fpsyt.2023.1188101/full#supplementary-material
References
[References listed in original article, numbers 1-154]
Glossary
Keywords: Z-drugs, benzodiazepine binding sites, sex differences, adverse events, pharmacovigilance, side effects, FDA adverse event reporting system, benzodiazepine
Citation: Koniuszewski F, Vogel FD, Dajić I, Seidel T, Kunze M, Willeit M and Ernst M (2023) Navigating the complex landscape of benzodiazepine- and Z-drug diversity: insights from comprehensive FDA adverse event reporting system analysis and beyond. Front. Psychiatry. 14:1188101. doi: 10.3389/fpsyt.2023.1188101
Received: 16 March 2023; Accepted: 05 June 2023; Published: 23 June 2023.
Edited by:
Yasser Khazaal, Université de Lausanne, Switzerland
Reviewed by:
Esa R. Korpi, University of Helsinki, Finland Hilary S. McCarren, US Army Medical Research Institute of Chemical Defense, United States
Copyright © 2023 Koniuszewski, Vogel, Dajić, Seidel, Kunze, Willeit and Ernst. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY).
