Investigation of [3H]diazepam derivatives as allosteric modulators of GABAA receptor α1β2γ2 subtypes: combination of molecular docking/dynamic simulations, pharmacokinetics/drug-likeness prediction, and QSAR analysis.

In this paper, a data set of [3H] diazepam derivatives was analyzed using various computational methods: molecular docking/dynamic simulations, and QSAR analysis. The main aims of these studies are to understand the binding mechanisms by which benzodiazepines allosterically modulate GABAA receptor α1β2γ2 subtypes, from inducing neuronal inhibition at lower doses to the anesthetic effect at higher doses, and also, to define the structural requirements that contribute to improving the response of GABAA/α1β2γ2 receptor to benzodiazepine drugs. The results of the molecular docking study allowed selecting Ro12-6377 and proflazepam as the best modulators for the four binding sites simultaneously. Subsequently, the stability of the selected complexes was investigated by performing molecular dynamics simulation. The latter confirmed the features of both modulators to exert direct effects on the chloride-channel lining residues. Pharmacokinetics and drug-likeness profile were assessed through in silico tool. Furthermore, a QSAR analysis was conducted using an improved vemolecular dynamics simulations proposed byrsion of PLS regression. The goodness of fit and the predictive power of the resulting PLS model were estimated according to internal and external validation parameters: R 2 = 0.632, R 2 adj = 0.584, F = 12.806; p-value = 6.2050e - 07, Q 2 loo = 0.639, and Q 2 F3 = 0.813. Clearly, the obtained results ensure the predictive ability of the developed QSAR model for the design of new high-potency benzodiazepine drugs. Supplementary Information: The online version contains supplementary material available at 10.1007/s11224-022-02029-4.

Introduction

Gamma-aminobutyric acid type A (GABAA) are fast-acting ionotropic receptors that belong to the superfamily of Cys loop-type ligand-gated ion channels (PLGICs) [1]. This type of GABA receptor is found at 20–50% of brain synapses [2].GABAARs assemble into pentameric isoforms consisting of five heteromeric subunits that surround a selective chloride-conducting pore. The general architecture of each mature subunit is constructed by the sequence of 450 amino acid residues [3]. In total, 200–250 amino acids contribute in the extracellular domain (ECD) to form a long hydrophilic N-terminal α-helix followed by ten β-strands folded into a β-sandwich containing the Cys-loops. A total of 85–255 amino acids contribute in the transmembrane domain (TMD) to form four membrane α-helices (named from TM1 to TM4) connected by three loops (short intracellular loop links TM1-TM2, short extracellular loop links TM2-TM3, and long cytoplasmic loop links TM3-TM4) and terminating with one small extracellular C-terminal [4]. The five TM2α helices were assembled in the center to form the chloride-channel liner, with a possible contribution fromTM1 [5]. Upon activation, GABAARs allow negative chloride ions to flow through the ion channel into the cell, thereby inhibiting neuronal excitability for a short time (phasic inhibition) or for a long time (tonic inhibition) [3, 6].

In the mammalian brain, GABAARs are assembled from 8 different families of subunits containing a total of 19 subtypes: α1-6, β1-3, γ1-3,δ, ε, π, θ, and ρ1-3. The predominant GABAARs isoform consists of two -subtypes, two β-subtypes, and one γ-subtype, or one δ-subtype. γ-containing receptors distribute mainly in the synaptic sites and account for approximately 90% of GABAARs in the adult brain. In contrast, δ-containing receptors are abundant in the extrasynaptic sites located in specific brain areas such as the hippocampus, amygdala, neocortex, thalamus, hypothalamus, and the cerebellum [6]. Owing to its complex subunit assembly, GABAARs contain large numbers of allosteric binding sites that make them the most important drug target in the central nervous system (CNS). They modulate by two endogenous molecules (neurosteroids and endocannabinoid 2-arachidonoylglycerol) and a wide range of exogenous molecules; the most important among them are benzodiazepines(BDZs), barbiturates, the intravenous general anesthetics propofol and etomidate, alcohols, the competitive antagonist bicuculline, and the channel blocker picrotoxin [3].

The α1β2γ2 combination constitutes approximately 43–60% of GABAARs in the adult brain [6], which makes it the subject of many previous reports attempting to provide high-resolution structural data that illuminate atomic mechanisms of drug recognition. The most important among them are the cryo-electron microscopy structures complemented with mutagenesis, electrophysiology, and molecular dynamics simulations proposed by Kim et al. [2, 7]. Kim et al. demonstrated, in addition to the classical BDZ site at the ECD / interface, the presence of three additional binding sites in the TMD: two at / interfaces and a third at the / interface (Fig. 1). / sites were also previously observed by Masiulis et al. [8] through their study on α1β3γ2Rs. The pharmacological role of BDZs exerts in the presence of two GABA in their binding sites at the ECD / interfaces. More precisely, BDZs react as positive allosteric modulators that enhance the affinity of GABARs for GABA, leading to an increase in the conduction of chloride ions through the ionic channel into the cell, hyperpolarization of neurons, and thus inhibition of neuronal excitability. The BDZ binding loci at both TMD interfaces were identified to be also the binding sites for etomidate and propofol. Also, the TMD binding locus located at the / interface was estimated to overlapped in part with that of the phenobarbital. Etomidate, propofol, and barbiturate are general intravenous anesthetics that differ from the benzodiazepines in their ability—at higher doses—to directly activate GABAARs without the need for GABA. In a mechanism similar to that of barbiturate, occupation of the / interface by diazepam (DZP) leads to close the gap presented at the interface between the two γ2-β2 subunits, which according to Kim et al., may explain the anesthetic properties observed during the administration of a high dose of DZP. Currently, the [3H]diazepam derivatives are the most common psychoactive drugs used to treat epilepsy, insomnia, muscle spasms, anxiety, alcohol withdrawal, and panic disorder [9].

Fig. 1

Binding sites of the two endogenous agonists (GABA) and diazepam. a ECD interface for the binding of diazepam (a) and GABA (a’ and a’’). b The three TMD interfaces identified for the binding of diazepam [7]

The strategy of combining molecular docking, molecular dynamics simulation, and QSAR analysis has emerged as a practical tool in the process of drug development through computational techniques. Its main advantage lies in improving the success rate of drug screening in less time and at a lower cost [10]. Molecular docking simulations are designed to determine the best ligand/target binding mode that generates the biological response. In practice, it allows the screening of large libraries of compounds by implementing fast and inexpensive docking algorithms. Its process is based on the production of all possible ligand poses within the binding target, and associates each pose with a score value that approximates its free energy landscape. At the ends of simulations, the best binding modes are ranked based on the values of the latter, and thus the most appropriate complexes are selected. Mostly, after molecular docking simulations, the best-docked complexes are subject to stability investigation through molecular dynamics simulation. Monitoring the dynamic profile of complexes over a certain time range provides the advantage of detecting various internal motions and conformational changes that occur in the binding site. Hence, validate the docking protocols. Otherwise, MD simulations can be generated before performing molecular docking for several objectives such as optimizing the target structure and ensuring its flexibility, quantifying the ligands/target free binding energies,…etc., as well as, during the docking process, to accurately detect the binding locus and properly dock ligands [11].

Quantitative structure–activity relationships (QSARs) are statistically derived models mainly devoted to predicting the activities of new chemical entities from knowledge of their chemical structures. QSAR models quantitatively correlate the physicochemical and biological properties of compounds with their biological responses [12]. The success of any QSAR model depends on many factors, including the selection of practical statistical techniques for model development and validation strategies. The quality of models is evaluated through an internal validation process usually based on the use of the cross-validation method (CV), whereas the predictive power can be estimated using independent test data that was not involved in model generation [13]. In the literature, considerable researchers have sought to study the BDZ by performing QSAR analysis. Against this background, fifty-seven compounds from the dataset of interest in our study were examined by Maddalena and Johnston [14] using the back-propagation artificial neural network method (ANN) and multilinear regression analysis (MLR). The two-layer ANN model gave an excellent correlation between the binding affinities of the 38 compounds of the training set (Rt = 0.941) and the 10 selected input variables (π7, F7, MR1, MR2', R1, F2, MR6, μ1, δp8, and δm3), where F, R, μ, δp, and δm denote polar constant, resonance constant, dipole moment, hammett-para constant, and hammett-meta constant, respectively. The predictive power was tested with an external set of 19 compounds, and an optimal cross-validation coefficient (Rcv = 0.910) was found. Later, So et al. [15] re-examined the data provided by Maddalena et al. with an improved version of the genetic neural network (GNN). The top-ranking variables selected by GNN shared only four variables with the ANN selection (π7, F7, MR1, and MR2').Therefore, an accurate comparison and a further discussion were made between the results of the two statistical methods, as well as the three highly predictive models (T6-2 # 1–3) which were combined with the optimal functional groups proposed by Maddalena to be placed in positions 1, 7, and 2' and used to design 20 new BDZ derivatives with their predicted activities.

This research aims to highlight the binding mechanism by which a data set of classical benzodiazepines allosterically modulates GABAA receptor α1β2γ2 subtypes. In addition to the well-known ECD binding interface, our study was further expanded to include the three TMD interfaces that were recently discovered. In the first step, a data set of [3H]diazepam derivatives was subjected to fast screening through the molecular docking approach. Subsequently, molecular dynamics simulation and pharmacokinetics/drug-likeness evaluations were performed to refine the best-docked complexes. Finally, an improved version of PLS regression was implemented to quantify structural features that contribute to the improvement of BDZ/α1β2γ2Rs responses. Throughout this paper, the results have been interpreted in light of the combination of the cited approaches.

Materials and methods

Biological data

The results in vitro for the 50% inhibition of the binding of [3H]diazepam to homogenates of rat brain cell membranes by BDZs expressed as log (1/C) reported earlier in the review of Hadjipavlou-Litina and Hansch [16] (Table 1) were investigated to perform a molecular docking simulation and to predict QSAR model using PLS analysis. According to Micheli et al. [17], the good structural diversity allows this dataset to be optimal for undergoing QSAR analysis.

Table 1

The classical BDZ data set under study [16]

N°	Name	R1	R3/ R5/R6	R7	R8/ R9	R2’	R6’	Log (1/C)obs
1	Ro05-4318/ Ro05-3418	CH3	R3,R5,R6 = H	NH2	R8,R9 = H	H	H	6.34
2	Ro05-3072	H	R3,R5,R6 = H	NH2	R8,R9 = H	H	H	6.41
3	Ro05-4528	CH3	R3,R5,R6 = H	CN	R8,R9 = H	H	H	6.42
4	Ro05-2921	H	R3,R5,R6 = H	H	R8,R9 = H	H	H	6.45
5	Ro20-7736	CH3	R3,R5,R6 = H	NHOH	R8,R9 = H	F	H	7.02
6	Ro05-4619	H	R3,R5,R6 = H	NH2	R8,R9 = H	Cl	H	7.12
7	Ro20-5397	H	R3,R5,R6 = H	CHO	R8,R9 = H	H	H	7.37
8	Ro05-3061	H	R3,R5,R6 = H	F	R8,R9 = H	H	H	7.40
9	Ro20-2533	H	R3,R5,R6 = H	C2H5	R8,R9 = H	H	H	7.44
10	Ro20-2541	CH3	R3,R5,R6 = H	CN	R8,R9 = H	F	H	7.52
11	Ro20-5747	H	R3,R5,R6 = H	CHCH2	R8,R9 = H	H	H	7.62
12	Ro05-4336	H	R3,R5,R6 = H	H	R8,R9 = H	F	H	7.68
13	Ro20-3053	H	R3,R5,R6 = H	COCH3	R8,R9 = H	F	H	7.74
14	Triflunordazepam	H	R3,R5,R6 = H	CF3	R8,R9 = H	H	H	7.89
15	Diazepam	CH3	R3,R5,R6 = H	Cl	R8,R9 = H	H	H	8.09
16	Ro07-5220	CH3	R3,R5,R6 = H	Cl	R8,R9 = H	Cl	Cl	8.26
17	Ro14-3074	H	R3,R5,R6 = H	N3	R8,R9 = H	F	H	8.27
18	Flunitrazepam	CH3	R3,R5,R6 = H	NO2	R8,R9 = H	F	H	8.42
19	Ro05-3590	H	R3,R5,R6 = H	NO2	R8,R9 = H	CF3	H	8.45
20	Norflurazepam	H	R3,R5,R6 = H	Cl	R8,R9 = H	F	H	8.7
21	Delorazepam	H	R3,R5,R6 = H	Cl	R8,R9 = H	Cl	H	8.74
22	Clonazepam	H	R3,R5,R6 = H	NO2	R8,R9 = H	Cl	H	8.74
23	Fonazepam	H	R3,R5,R6 = H	NO2	R8,R9 = H	F	H	8.82
24	Ro05-6822	CH3	R3,R5,R6 = H	F	R8,R9 = H	F	H	8.29
25	Ro05-4865	CH3	R3,R5,R6 = H	F	R8,R9 = H	H	H	7.77
26	Ro05-6820	H	R3,R5,R6 = H	F	R8,R9 = H	F	H	8.13
27	Nordazepam	H	R3,R5,R6 = H	Cl	R8,R9 = H	H	H	8.03
28	Ro07-3953	H	R3,R5,R6 = H	Cl	R8,R9 = H	F	F	8.79
29	Difludiazepam	CH3	R3,R5,R6 = H	Cl	R8,R9 = H	F	F	8.39
30	Ro07-5193	H	R3,R5,R6 = H	Cl	R8,R9 = H	Cl	F	8.52
31	Ro22-3294	H	R3,R5,R6 = H	Cl	R8,R9 = H	Cl	Cl	8.15
32	Nitrazepam	H	R3,R5,R6 = H	NO2	R8,R9 = H	H	H	7.99
33	Methylclonazepam	CH3	R3,R5,R6 = H	NO2	R8,R9 = H	Cl	H	8.66
34	7-Aminoflunitrazepam	CH3	R3,R5,R6 = H	NH2	R8,R9 = H	F	H	7.19
35	Ro12-6377	CH3	R3,R5,R6 = H	NHCONHCH3	R8,R9 = H	F	H	6.34
36	Halazepam	CH2CF3	R3,R5,R6 = H	Cl	R8,R9 = H	H	H	7.04
37	Pinazepam	CH2CCH	R3,R5,R6 = H	Cl	R8,R9 = H	H	H	7.03
38	Prazepam	CH2C3H5	R3,R5,R6 = H	Cl	R8,R9 = H	H	H	6.96
39	Motrazepam	CH2OCH3	R3,R5,R6 = H	NO2	R8,R9 = H	H	H	6.37
40	Ro20-1310	C(CH3)3	R3,R5,R6 = H	Cl	R8,R9 = H	H	H	6.21
41	Ro07-2750	CH2CH2OH	R3,R5,R6 = H	Cl	R8,R9 = H	F	H	7.61
42	Ro08-9013	(CH2)2OCH2CONH2	R3,R5,R6 = H	Cl	R8,R9 = H	F	H	7.37
43	Proflazepam	CH2CHOHCH2OH	R3,R5,R6 = H	Cl	R8,R9 = H	F	H	6.85
44	Ro22-4683	C(CH3)3	R3,R5,R6 = H	NO2	R8,R9 = H	Cl	H	6.52
45	Ro11-4878	H	R3 = (s)CH3 R5,R6 = H	Cl	R8,R9 = H	F	H	8.46
46	Meclonazepam	H	R3 = (s)CH3 R5,R6 = H	NO2	R8,R9 = H	Cl	H	8.92
47	Ro11-6896	CH3	R3 = (s)CH3 R5,R6 = H	NO2	R8,R9 = H	F	H	8.15
48	L48	CH3	R3 = (rac)CH3 R5,R6 = H	Cl	R8,R9 = H	H	H	7.31
49	Temazepam	CH3	R3 = (rac)OH R5,R6 = H	Cl	R8,R9 = H	H	H	7.79
50	L50	CH3	R3 = (rac)Cl R5,R6 = H	Cl	R8,R9 = H	F	H	8.27
51	L51	H	R3,R5 = H R6 = CH3	CH3	R8,R9 = H	H	H	6.77
52	Ro07-4419	H	R3,R5,R6 = H	H	R8,R9 = H	F	F	7.72
53	Ro05-4520	CH3	R3,R5,R6 = H	H	R8,R9 = H	F	H	7.85
54	Ro05-4608	CH3	R3,R5,R6 = H	H	R8,R9 = H	Cl	H	8.42
55	Ro05-3546	H	R3,R5 = H R6 = Cl	H	R8,R9 = H	H	H	6.49
56	Ro13-0699	CH3	R3,R5 = H R6 = Cl	H	R8,R9 = H	F	H	6.82
57	Ro07-6198	H	R3,R5,R6 = H	H	R8 = Cl R9 = H	F	F	7.55
58	Ro20-8895	H	R3,R5,R6 = H	H	R8 = CH3 R9 = H	F	H	7.72
59	Ro13-0593	CH3	R3,R5,R6 = H	H	R8 = H R9 = Cl	F	H	7.14
60	L60	CH3	R3,R5 = H R6 = Cl	H	R8 = Cl R9 = H	F	H	6.52
61	Ro22-6762	CH3	R3,R5,R6 = H	Cl	R8 = Cl R9 = H	H	H	7.40
62	Ro20-8065	H	R3,R5,R6 = H	Cl	R8 = Cl R9 = H	F	H	8.44
63	Ro20-8552	H	R3,R5,R6 = H	CH3	R8 = Cl R9 = H	F	H	7.85
64	L64	H	R3,R5,R6 = H	Cl	R8 = H R9 = Cl	H	H	7.43
65	L65	H	R3,R5,R6 = H	Cl	R8 = H R9 = CH3	H	H	7.28

Molecular descriptor generation

First, a gradient norm limit of 0.1 kcal/(Å mol) was chosen for the pre-optimization of sixty-five BDZ derivatives, using the molecular mechanics force field (MM+) method included in HyperChem package version 8.08 [18]. Then, the molecular geometries were optimized at the DFT/Ub3lyp/6–311++G(d,p) level of theory using the Gaussian 09 W software [19]. For all stationary points, there is no imaginary frequency at the optimized molecular geometries ensuring that the optimized structures are at the minimum on the potential energy surface. The atomic charges qN1, qC3, qN4, qC6, qC7, qC8, qC9, qC2', qC6', and the dipole moment (DM) have been assessed using the ChelpG electronic population scheme [20]. The number of hydrogen-donors (HD), number of hydrogen-acceptors (HA), molecular lipophilicity (Log P), and molar refractivity (MR) have been computed using MarvinSketchversion 20.21 [21]. The flexible torsions (FT) have been computed using Molegro Virtual Docker version 5.5 [22]. Finally, the hydrophobic constants in positions 7 (πC7) and 2’ (πC2') have been extracted from the literature [23].

QSAR analysis

PLS regression

The chemical, physical, topological, and quantum properties are necessarily correlated for a given molecule. This is a reflection of the innate properties of the system, and additional data collected in the same way will show the same collinearity [24]. Indeed, PLS regression is a useful method for multivariate data containing correlated molecular descriptors. This method based on dimension reduction technique builds orthogonal components, often called factors or latent variables, as linear combinations of the original predictor variables [25]. PLS constructs these components while considering the observed response values, leading to a parsimonious model with reliable predictive power [26].

In this work, the number of components used in PLS is chosen by a fivefold cross-validation method [27]. The dataset is randomly divided into training dataset (80%) and testing dataset (20%). Training sets were used for model development and test sets for model external validation. Before conducting PLS analysis, each response variable is scaled to unit variance by dividing it by its standard deviation, and the molecular descriptors are centered by subtracting the average value and scaled to unit variance. First, variable selection by stepwise regression method is used to identify the best subset of molecular descriptors. This is a combination of backward and forward selection [25]. The objective is to use the minimum number of descriptors to develop a good predictive model. Thus, we must select the good subsets of descriptors. However, it should be noted that the subset of molecular predictors that do the best at meeting well-defined objective criteria can be highly variables depending on precisely which observations are included in the training set. In addition, the best training model does not necessarily guarantee a better quality of prediction. This depends on the training and test sets obtained from the original dataset. For this reason, we conducted a statistical simulation for which 10,000 splits were performed resulting in 10,000 training and test sets. For each simulation, regression diagnostics for detecting possible outliers was carried out by computing leverage values (h) for identifying outlying x-variables and studentized deleted residuals for identifying outlying -variables. Note that once outliers have been detected the model is regenerated excluding the outlying observations from the dataset. Thereafter, the best model is selected on each training set resulting in 10,000 best training models following the Bayesian Information Criterion (BIC) [28]. This criterion is chosen because it penalizes larger models more heavily and will tend to select a smaller subset of descriptors in comparison to other criteria [29]. The best choice of descriptors will balance fit with model size. Subsequently, among these 10,000 best models, we sought to select the best molecular descriptors according to the highest probability of their occurrence.

To verify a model’s predictive ability, the developed QSAR model is quantified using the coefficient of determination (R2) [30], the adjusted coefficient of determination [31], and the Fisher statistics (F) [32]. This latter is computed to judge the overall significance of the regression model. The external predictive ability of the developed QSAR model is determined by computing the leave-one-out cross-validation coefficient [32] and the predictive squared correlation coefficient [33-35]. The external validation ensures the predictability of the developed QSAR model for the prediction of untested molecules [13].

Molecular docking protocol

The electron microscopy structure of the human GABAA receptor α1β2γ2 subtypes in complex with GABA plus the DZP structures (PDB ID:6X3X, Resolution = 2.92 Å) was downloaded from RCSB Database (http://www.rcsb.org). The downloaded PDB file contains nine chains: five chains denote the subunits α1(B and D), β2(A and C), and γ2 (E), and four chains denote the Fab-chains (named from I to K).

MOE 2014.0901software package [36] was used to prepare the four benzodiazepine binding sites: the classical site at the ECDD+/E− chain interface and the three TMD sites at the A+/B−, C+/D−, and E+/A− chain interfaces. For each binding site, all co-crystallized ligands and non-essential subunits were removed from the α1β2γ2-DZP complex. Then, after structure correction, protonation at neutral medium (PH = 7), and cavity detection, the native co-crystallized DZP structure was re-docked into the selected binding site pocket. The method was validated by giving the best binding pose which has a low RMSD value (root-mean-square deviation).

The BDZ structures previously optimized using the DFT method were converted into databaseinput files and docked one by one into the four DZP binding pockets using the semi-flexible docking process of MOE 2014.0901 package [36]. During the process, the conformation of the receptor was fixed, while ligands remained flexible. Here, the best binding poses were selected according to the lowest energy score values, registered in the PDB file, and visualized using BIOVIA Discovery Studio visualizer v20.1.0.19295 package [37].

Molecular dynamics protocol

The best-docked ligand/α1β2γ2 complexes were subjected to stability tests using molecular dynamics simulation (MD). MD simulation was implemented using the “compute” option included in the MOE 2014.0901 software [36]. First, the selected complexes were prepared by deleting the DZPco-crystallized structure, fixing hydrogens, and fixing charges. Thereafter, the parameters of the “dynamics” tool were adjusted to execute the combination Nosé-Poincaré-Andersen (NPA) algorithm/Merck molecular force field (MMFF94x) [38, 39], with enabling bonding, van der Waals, electrostatics, and restraints. The protocols were settled for an equilibrium period of 100 ps followed by a production period of 900 ps, at a constant temperature of 310 K. Finally, the variations in potential energies, U (kcal/mol), as a function of time, t (ps), are retained and plotted using Origin 6.0 software [40].

Pharmacokinetics/drug-likeness prediction

In silico estimation of pharmacokinetic properties and prediction of drug-likeness were carried out by using the free web tool SwissADME [41]. Our study is based on the prediction of the following pharmacokinetic parameters: gastrointestinal absorption (GI), P-glycoprotein (P-gp) substrate, blood–brain barrier (BBB) penetration, and cytochrome enzyme (CYP) inhibition. Indeed, out of 57 human CYP450 enzymes, the CYP1A2, CYP2C9, CYP2C19, CYP2D6, CYP3A4, and CYP2E1 metabolize 90% of drugs [42]. In addition, the drug-likeness prediction is based on several rules such as Lipinski, Ghose, Veber, Egan, and Muegge.

Results and discussion

Molecular docking simulation

Orientations, interactions, and binding affinities of 65 positive allosteric modulators of GABAARs (BDZ dataset, Table 1) were investigated in four distinct BDZ binding sites: the classical site at the ECD / interface, the TMD sites at the two / and / interfaces and the TMD site at the / interface (Fig. 1). Residues involved in each active pocket were detected using the site finder wizard implemented in the MOE 2014.0901 software [36] (Table 2). During the docking process, the conformation of residues remains unchanged, while ligands are altered. The best binding modes of re-docked DZPs were selected based on the given root-mean-square deviation values (RMSD). Generally, docking protocols that are able to generate the same co-crystallized binding modes with an RMSD value of less than 1.5 or 2 A° (depending on ligand size) or even better, less than 1 A° are considered validated [43, 44], whereas the best binding modes of docked ligands were selected according to the given Sscore values. The binding free energy score (Sscore) is a quantitative estimate of the most stable binding pose between the target macromolecule and ligand. Among the generated poses, the best ones are those with the most negative energy score values.

Table 2

Residues involved in each active pocket

Binding site		Residues involved in active pockets
ECD	α1+/γ2- interface (classical binding site)	1: (PHE100 PHE101 HIS102 ASN103 GLU138 PRO140 PRO154 LYS156 SER159 TYR160 ALA161 VAL203 GLN204 SER205 SER206 THR207 TYR210) 2: (ASP56 MET57 TYR58 ASN60 SER61 ASP75 PHE77 ALA79 MET130 THR142 ARG144 SER186 GLU189 ASP192 SER195)
TMD	β2+(A)/α1-(B) interface	1: (ILE255 VAL258 LEU259 MET261 THR262 ASN265 THR266 ARG269 GLU270 ASP282 LEU285 MET286 PHE289 VAL290) 2: (ILE228 GLN229 LEU232 PRO233 MET236 THR237 LEU240 PHE258 THR261 THR262 THR265 LEU269 SER272)
	β2+(C)/α1-(D) interface	1: (MET261 THR262 ASN265 THR266 ARG269 ASP282 LEU285 MET286 PHE289 VAL290) 2: (VAL227 ILE228 LEU232 PRO233 MET236 THR237 THR265 LEU269)
	γ2+/β2- interface	1: (TYR220 PHE221 LEU223 GLN224 THR225 MET227 PRO228 LEU231 ILE232 THR263 ILE264 HIS267 LEU268 THR271 LEU272) 2: (MET276 THR277 SER280 THR281 ALA283 ARG284 LYS285 LYS289 MET296 ASP297 VAL300 SER301 PHE304 ILE305)

/

The energy score values of the 65 ligands docked in the /, /, /, and / interfaces are between (− 8.013 and − 6.046), (− 7.409 and − 5.566), (− 7.455 and − 5.463), and (− 7.546 and − 5.425), respectively (Table S1, Supplementary materials). As it evident, the studied ligands have rather scattered affinities around the reference values: − 7.003, − 6.159, − 6.261, and − 6.347, respectively.

Surveying the first five ligands having the lowest energy scores (Table 3), the lowest Sscore value in the three /, /, and / sites was assigned to Ro12-6377. Also, Ro12-6377 exhibits the second-lowest Sscore value in the / interface. Therefore, out of the 65 studied ligands, Ro12-6377 is the ligand that exhibited the highest predicted affinity towards the four sites, simultaneously. By similar reasoning, the second highest predicted affinity was rated to proflazepam.

Table 3

The first five ligands having the highest binding affinity for the four binding interfaces

ECD α1+/γ2- interface		TMD β2+(A)/α1-(B)interface		TMD β2+(C)/α1-(D) interface		TMD γ2+/β2- interface
Ligand	Sscore	Ligand	Sscore	Ligand	Sscore	Ligand	Sscore
Ro12-6377	− 8.013	Ro12-6377	− 7.409	Ro08-9013	− 7.455	Ro12-6377	− 7.546
Meclonazepam	− 7.856	Proflazepam	− 7.170	Ro12-6377	− 7.265	Ro08-9013	− 7.070
Methylclonazepam	− 7.831	Ro07-2750	− 7.152	Proflazepam	− 7.077	Pinazepam	− 6.990
Proflazepam	− 7.807	Pinazepam	− 7.146	Pinazepam	− 6.848	Proflazepam	− 6.972
Ro11-6896	− 7.785	Motrazepam	− 7.087	Meclonazepam	− 6.774	Ro07-2750	− 6.918

By comparing the affinities of Ro12-6377 towards the four sites, the / site was defined to be the principal target for Ro12-6377. So, when a high dose of Ro12-6377 is administered, it mainly acts to bind at the ECD / interface as it is the high-affinity binding site, and then, acts to bind at the TMD / and the two / interfaces as they are the second and third high-affinity binding sites, respectively. On the other hand, proflazepam tends to fill thetwo / interfaces before moving to occupy the / site. Recently, from a database of 7922 compounds, proflazepam was selected among the top-100 docked ligands at the binding pocket of SARS-CoV-2 Main Protease crystallized in holo-form [45].

Likewise, the affinities of the remaining 63 ligands (Table S1, Supplementary materials) towards the four binding sites were compared. Unexpectedly, our results are inconsistent with previous findings indicating that the classical site is always the main target of all classical benzodiazepines. As we can see, / and / are expected to be, respectively, the main binding sites for (Ro05-4336, Ro07-2750, Ro05-3546) and (L51,Ro13-0699). Moreover, in some cases, both ECD and TMD sites share the same binding affinity towards the bound ligand, as is evident for Ro05-2921, Ro20-5397, Pinazepam, Ro07-4419, Ro05-3546, and Ro20-8552, which makes it difficult to distinguish, accurately, the main target for the binding. Consequently, after the discovery of the three TMD sites, it became necessary to expand previous findings indicating that the main target of classical benzodiazepines is always located at the ECD / interface.

The binding modes of both ligands in each of the four binding interfaces are located at the same level as the co-crystallized DZP (Figs. S1 and S2, Supplementary materials). However, their docking orientations are markedly not equivalents, with the exception of the binding orientation of proflazepam in the / interface where it shares an almost perfect superimposition to the DZP-bound structure (Fig. S2 (d) in Supplementary materials). The favorable drug binding pose is determined by the distribution of polar and non-polar regions along the surface of the ligand and its complementary target binding site. Thus, while the non-polar regions create hydrophobic interactions mainly contributing to the binding affinity of the drug towards the biological target, polar regions create electrostatic points contributing to modulating the drug binding kinetics (specificity and orientation) [46]. In order to estimate all possible interactions, the docking-outputs generated by MOE software were converted into (.pdb) files and visualized with the default parameters of BIOVIA DS visualizer v20.1.0.19295 package [37]. As can be seen, the binding interactions for both ligands with the four target-sites residues exhibit the formation of four types of interactions, most of which are of type hydrogen bonds and hydrophobic interactions. The standard values of distances and energy cutoffs for considering the formation of hydrogen bonds with specific target residues are categorized into three subtypes: strong bonds (2.2–2.5 Å, E: 14–40 kcal/mol), moderate bonds (2.5–3.2 Å, E: 4–15 kcal/mol), and weak bonds (3.2–4.0 Å, E > 4 kcal/mol) [47]. Generally, strong bonds in targeted-ligand interactions are undesirables as they mostly tend to have a covalent character. This later hinders the process of drug liberation from its receptor, thus increasing the risk of drug toxicity. Nevertheless, new insights on how to address the pharmacological advantages/potential risks balance of covalent drugs have been emerged and discussed in the literature [48, 49]. Referencing to this latter, several standard values of distance cutoffs to consider the formations of hydrophobic interactions have been found. Janiak [50] suggested that the optimum range is between 3.3 and 3.8 A°, while other researchers have suggested a relatively higher range [51-53].

Owing to the lack of previous experimental and theoretical studies on Ro12-6377 and proflazepam, we will attempt to suggest mechanisms for how they modulate α1β2γ2Rs signaling by binding at the ECD and the TMD interfaces, based on the results obtained from molecular docking and molecular dynamics simulation.

Allosteric modulation of the classical binding site

The bulky structure of Ro12-6377, compared toDZP, enabled the NHCONHCH3 group attached at C7 and both phenyls (B) and (C) to penetrate deeper into the binding site, whereas the diazepine ring exerts its effect by interacting with the residues located in front of the pocket (Figs. 2a and 3a and Table 4). The methyl group attached at N1 forms two hydrophobic interactions type Alkyl-Alkyl and Pi-Alkyl with α1Val203 andγ2Tyr58, respectively. Concurrently, the π-electron clouds of γ2Tyr58 and γ2Phe77 were involved in two Pi-Pi stacked interactions with the π-electron cloud of phenyl (B). The linear backbone skeleton of NHCONHCH3 group formed an Alkyl-Alkyl interaction with the side chain of γ2Ala79, received two moderate H-bonds from α1Ser206side chain, and gave two moderate H-bonds to the main carbonyl of γ2Phe77. The four moderate H-bonds lead to the formation of two intermolecular pentameric rings that contribute to the stability of Ro12-6377 at the binding site. Also, in addition to the strongH-donor bond received from α1His102, an intramolecular interaction was observed between the fluorine atom bound at C2' andC2. This interaction could be explained by the strong withdrawing property of the oxygen that madeC2 a more electrophilic center able to receive the nucleophilic attacks. The pendant phenyl (C) is involved deeper into the DZP pocket where it is delineated by the hydrophobic side chains of γ2Phe77, α1Phe100, α1His102, α1Ter160, and α1Ter210. Here, a hydrophobic interaction type Pi-Pi T-shaped was observed between the π-electron cloud of ring (C) and the surrounded side chain of α1Phe100.

Fig. 2

Binding modes resulting from the molecular docking of Ro12-6377and proflazepam at the interfaces of a ECD /, b TMD /, c TMD /, and d TMD /

Fig. 3

Binding interactions resulting from the molecular docking of Ro12-6377 at the interfaces of a ECD /, b TMD /, c TMD /, and d TMD /

Table 4

Detailed binding interactions resulting from the molecular docking of co-crystallized DZP, Ro12-6377, and proflazepam (PLZ) at the classical site

ECD α1+/γ2- interface
Ligand name	Posita	Category	Type of interactions	From	From chemistry	To	To chemistry	Distb (A°)
DZP	2	Hydrogen bond	Conventional H-bond	D:SER206:H	H-donor	D:DZP404:O1	H-acceptor	2.68
	4			D:DZP404:H151		D:SER205:OG		2.55
	2		C-H bond	D:SER205:HA		D:DZP404:O1		2.65
	1			D:DZP404:H202		D:GLN204:O		2.94
	4	Hydrogen bond; electrostatic	Pi-Cation; Pi-Donor H-bond	D:DZP404:H152	Positive; H-donor	E:PHE77	Pi-orbitals; Pi-orbitals	3.03
	Ring C	Hydrophobic	Pi-Pi stacked	D:TYR160	Pi-orbitals	D:DZP404	Pi-orbitals	5.24
				D:TYR210		D:DZP404		4.22
				E:PHE77		D:DZP404		5.49
			Pi-Pi T-shaped	D:PHE100		D:DZP404		5.36
	1		Alkyl	D:DZP404:C20	Alkyl	D:VAL203	Alkyl	5.21
	7		Pi-Alkyl	D:HIS102	Pi-orbitals	D:DZP404:CL		4.32
	1			E:TYR58		D:DZP404:C20		4.37
	7			E:PHE77		D:DZP404:CL		5.12
	Ring B			D:DZP404		D:VAL203		5.24
Ro12-6377	7	Hydrogen bond	Conventional H-bond	D:SER206:HG	H-donor	:*0:O	H-acceptor	2.57
				:*0:H		E:PHE77:O		3.06
	2’		C-H bond	D:HIS102:HE1		:*0:F		2.47
	7			D:SER206:HB3		:*0:O		2.70
				:*0:H		E:PHE77:O		2.98
	Internalc	Halogen	Halogen (Fluorine)	:*0:C	Halogen acceptor	:*0:F	Halogen	3.66
	Ring B	Hydrophobic	Pi-Pi stacked	E:TYR58	Pi-orbitals	:*0	Pi-orbitals	3.93
				E:PHE77		:*0		5.00
	Ring C		Pi-Pi T-shaped	D:PHE100		:*0		5.23
	7		Alkyl	E:ALA79	Alkyl	:*0:C	Alkyl	3.55
	1			:*0:C		D:VAL203		4.55
			Pi-Alkyl	E:TYR58	Pi-orbitals	:*0:C		4.88
PLZ	1	Hydrogen bond	Conventional H-bond	D:THR207:HG1	H-donor	:*0:O	H-acceptor	2.10
				D:TYR210:HH		:*0:O		2.85
	4		C-H bond	D:HIS102:HE1		:*0:N		2.57
	Ring B	Hydrophobic	Pi-Pi stacked	E:TYR58	Pi-orbitals	:*0	Pi-orbitals	4.73
				E:PHE77		:*0		4.55
	7		Pi-Alkyl	E:TYR58		:*0:CL	Alkyl	3.43

*0: ligand structure

aPosition of interaction

bdistance of interaction

cinternal interaction

Compared with Ro12-6377 and DZP, proflazepam penetrates less within the binding site (Figs. 2a and 4a and Table 4). Phenyl (C) is oriented outside the binding pocket, which explains the absence of any interaction on this pendant ring. Rings (A) and (B) seem to selectively accept rather than accept and donate bonds to the residues of subunits. The subunit acts through three H-bonds; two are oriented from the hydroxyl groups at the side chains of Thr207 and Tyr210 towards the first hydroxyl group of the CH2CHOHCH2OH substitute, and one is oriented from the side chain of His102 towards N4. The γ2 subunit prefers to act hydrophobically through two types of interactions: Pi-Pi stacked and Pi-Alkyl. The Pi-Pi stacked is created between the π-electron cloud of phenyl (B) and the Tyr58 and Phe77 side chains, while the Pi-Alkyl interaction links the chlorine atom at C7 to the side chain of Tyr58.

Fig. 4

Binding interactions resulting from the molecular docking of proflazepam atthe interfaces of a ECD /, b TMD /, c TMD /, and d TMD /

The binding modes of Ro12-6377 and proflazepam elucidated here agree with the previous finding indicating the importance of His102 in the recognition of classical benzodiazepines. With the exception of the His102Cys mutant in the α5 subunit, GABAARs that contain αHis102 mutation to any other residue suffer from total insensitivity to the DZP and its analogs [2]. The topological organization of the α4βγ2 and α6βγ2 receptors reveals the presence of a natural substitution of His102 by the Arg residue, which leads to steric problems affecting the binding of classical BDZ at their correspondent binding locus. This could explain the selectivity of classical BDZ towards αβγ2Rs containing the α1, α2,α3, and α5 subunits rather than those containing the α4 and α6 subunits [8, 9]. Else, the γ2Phe77Tyr mutant affects less the binding affinity of DZP but more strongly reduces that of its analogs containing the chlorine substitutes at the pendant phenyl (C). This finding was explained by the difference in flexibility between the two pendant phenyls since the presence of the chlorine atoms possibly caused unfavorable steric clashes with the side chain of the tyrosine residue. Furthermore, several other mutation findings were surveyed in detail in the previous researches [54-56].

Allosteric modulation of the three TMD binding sites

Despite the common topological organization between the combinations β2(A)/α1(B) and β2(C)/α1(D), the two TMD orthosteric pockets inserted at and interfaces are not qualitatively identical and thus may not be functionally equivalent. The binding pocket at interface was rated as the largest, and estimated to have a higher affinity for Ro12-6377 and proflazepam than the binding pocket at interface.

The binding modes of Ro12-6377 for both TMD interfaces are equivalents: the NHCONHCH3 groups and the fused benzodiazepine rings (A and B) are deeply embedded in the binding pockets, whereas the pendant phenyl rings (C) interact with the residues in front of the pockets (Fig. 2b, c). By superposing the binding modes, a rotation of 98.8° was observed between the two poses of the phenyl rings (C) (Fig. S3, Supplementary materials). The binding interactions at interface exhibit the formation of halogen interaction between the bound fluorine atom at C2' and N1 (Table 6). This intramolecular interaction leads the fluorine to orient towards the principal chain of β2Met286 where it stabilizes through the formation of a strong H-bond interaction (Fig. 3c), and consequently, leads the phenyl ring (C) to deviate by an angle of 98.8° from the phenyl plane observed at the interface (Fig. 3b). At this time, the binding interactions exhibit the formation of additional electrostatic and hydrophobic interactions between the substitute groups at N1 and C7 and the side chains of four residues from the β2 subunit: Val258, Leu285, Met286, and Phe289 (Table 5).

Table 6

Detailed binding interactions resulting from the molecular docking of co-crystallized DZP, Ro12-6377, and proflazepam (PLZ) at the TMD / interface

TMD β2+(C)/α1-(D) interface
Ligand Name	Posita	Category	Type of interactions	From	From chemistry	To	To chemistry	Distb (A°)
DZP	1	Hydrogen Bond	C-H bond	C:DZP406:H201	H-donor	D:ILE228:O	H-acceptor	3.08
				C:DZP406:H202		D:ILE228:O		2.56
	4	Hydrogen Bond; Electrostatic	Pi-Cation; Pi-Donor H-Bond	C:DZP406:H152	Positive; H-donor	C:PHE289	Pi-orbitals; Pi-orbitals	2.65
	Ring C	Hydrophobic	Pi-Sigma	D:PRO233:HB2	C-H	C:DZP406	Pi-orbitals	2.60
	1		Alkyl	C:DZP406:C20	Alkyl	C:MET286	Alkyl	4.68
				C:DZP406:C20		D:ILE228		4.76
				C:DZP406:C20		D:LEU232		4.82
				C:DZP406:C20		D:PRO233		4.46
	7			C:DZP406:CL		C:MET261		4.61
				C:DZP406:CL		C:LEU285		4.41
			Pi-Alkyl	C:PHE289	Pi-orbitals	C:DZP406:CL		4.99
	Ring C			C:DZP406		D:LEU269		5.50
	Ring B			C:DZP406		C:MET286		5.39
Ro12-6377	2’	Hydrogen bond; halogen	C-H bond; halogen (Fluorine)	C:MET286:HA	H-donor; halogen acceptor	:*0:F	H-acceptor; halogen	2.19
	3	Hydrogen bond	C-H bond	:*0:H	H-donor	D:ILE228:O	H-acceptor	2.95
				:*0:H		D:ILE228:O		3.03
	7			:*0:H		D:THR237:OG1		2.79
	Internalc	Halogen	Halogen (Fluorine)	:*0:N	Halogen acceptor	:*0:F	Halogen	3.50
	Ring B	Hydrophobic	Pi-Pi stacked	C:PHE289	Pi-orbitals	:*0	Pi-orbitals	4.14
	7		Alkyl	:*0:C	Alkyl	D:LEU240	Alkyl	4.63
	Ring B		Pi-Alkyl	:*0	Pi-orbitals	D:PRO233		4.40
	Ring C			:*0		C:MET286		4.54
				:*0		D:LEU232		4.94
				:*0		D:MET236		4.87
PLZ	2’	Hydrogen bond; halogen	Conventional H-bond; halogen (Fluorine)	C:ASN265:HD21	H-donor; halogen acceptor	:*0:F	H-acceptor; halogen	2.46
	4	Hydrogen bond	Conventional H-bond	C:ASN265:HD22	H-donor	:*0:N	H-acceptor	1.90
	1			:*0:H		D:ILE228:O		2.25
			C-H bond	D:PRO233:HD3		:*0:O		2.58
	Internalc			:*0:H		:*0:F		2.53
	Internalc			:*0:H		:*0:O		3.07
	2’	Halogen	Halogen (Fluorine)	D:ILE228:O	Halogen acceptor	:*0:F	Halogen	3.38
	Ring B	Hydrophobic	Pi-Pi stacked	C:PHE289	Pi-orbitals	:*0	Pi-orbitals	4.34
			Pi-Alkyl	:*0		D:PRO233	Alkyl	4.46
	Ring C			:*0		D:PRO233		4.52
				:*0		D:LEU269		5.19

*0: ligand structure

aPosition of interaction

bdistance of interaction

cinternal interaction

Table 5

Detailed binding interactions resulting from the molecular docking of co-crystallized DZP, Ro12-6377, and proflazepam (PLZ) at the TMD / interface

TMD β2+(A)/α1-(B) interface
Ligand Name	Posit a	Category	Type of interactions	From	From chemistry	To	To chemistry		Distb (A°)
DZP	1	Hydrogen bond	C-H bond	A:DZP406:H201	H-donor	B:ILE228:O	H-acceptor		2.93
				A:DZP406:H202		B:ILE228:O			2.89
	4	Electrostatic	Pi-Cation	A:DZP406:N15	Positive	A:PHE289	Pi-orbitals		4.22
	1	Hydrophobic	Alkyl	A:DZP406:C20	Alkyl	A:MET286	Alkyl		4.50
				A:DZP406:C20		B:ILE228			4.46
	7			A:DZP406:CL		A:MET261			4.42
				A:DZP406:CL		A:LEU285			4.79
			Pi-Alkyl	A:PHE289	Pi-orbitals	A:DZP406:CL			4.86
	Ring C			A:DZP406		B:PRO233			4.28
Ro12-6377	3	Hydrogen bond	C-H bond	:*0:H	H-donor	B:ILE228:O		H-acceptor	2.78
				:*0:H		B:ILE228:O			2.89
	1			:*0:H		A:LEU285:O			3.04
	7			:*0:H		B:THR237:OG1			2.51
	Ring B	Hydrophobic	Pi-Pi stacked	A:PHE289	Pi-orbitals	:*0		Pi-orbitals	4.11
	1		Alkyl	:*0:C	Alkyl	A:LEU285		Alkyl	4.64
				:*0:C		A:MET286			5.39
	7			:*0:C		A:VAL258			5.37
				:*0:C		B:LEU240			5.04
			Pi-Alkyl	A:PHE289	Pi-orbitals	:*0:C			5.30
	Ring B			:*0		B:PRO233			4.65
	Ring C			:*0		A:MET286			4.67
				:*0		B:LEU232			5.20
				:*0		B:MET236			4.66
PLZ	1	Hydrogen bond	Conventional H-bond	A:ARG269:HH12	H-donor	:*0:O		H-acceptor	2.10
				:*0:H		A:ASN265:OD1			1.97
	Internalc			:*0:H		:*0:O			2.44
	1		C-H bond	A:ARG269:HD3		:*0:O			2.72
	3			:*0:H		A:LEU285:O			2.90
	7	Hydrophobic	Alkyl	:*0:CL	Alkyl	B:LEU232		Alkyl	4.20
				:*0:CL		B:MET236			3.73
	Ring B		Pi-Alkyl	:*0	Pi-orbitals	A:MET286			4.91
				:*0		B:LEU232			5.42
	Ring C			:*0		B:PRO233			4.28

*0: ligand structure

aPosition of interaction

bdistance of interaction

cinternal interaction

Unlike Ro12-6377, the binding modes of proflazepam for both TMD and interfaces reveal significant differences in the docking orientations and binding interactions (Fig. 2b, c). At interface (Fig. 4b and Table 5), the chlorobenzene ring (B) was oriented to interact hydrophobically with side chains of residues located in front of the pocket:β2Met286, α1Met236, and α1Leu232. The pendant phenyl ring (C) was embedded deeper into the binding pocket, in a manner similar to that observed inDZP, and shares Pi-Alkyl interaction at a distance of 4.28A° with side chain of α1Pro233. The diazepine ring (A) established a moderate H-bond between one of the hydrogens of C3 and the main carbonyl of β2Leu285. The CH2CHOHCH2OH group was oriented towards the ECD where its first hydroxyl group was stabilized by strong intramolecular interaction type H-donor with lone pairs of the oxygen bond at C2, and with three intermolecular H-bonds formed with side chains of Arg269 and Asn265 of subunitβ2. At interface (Fig. 4c and Table 6), the pendant phenyl (C) is located at a higher level than rings (A) and (B), with its fluorine atom pointing towards the front of the binding site. This orientation drives the π-electron cloud into Pi-Alkyl interactions with side chains of Pro233 and Leu269 of subunit α1, as well as driven the Lone pairs of the fluorine atom to form halogen-halogen interaction with the main oxygen of α1Ile228 and receive a strong H-bond from the side chain of β2Asn265. This latter also shares a strong H-bond with N4. Phenyl (B) is located between the side rings of β2Phe289 and α1Pro233, which contributes to the formation of two hydrophobic interaction types Pi-Pi stacked and Pi-Alkyl. Its chlorine atom is oriented towards the bottom of the binding pocket, more precisely towards Thr265 of α1subunit. Here, no interactions were observed. The CH2CHOHCH2OH group exerts its influence by occupying the front of the binding pocket and sharing two H-bonds with α1Ile228 and α1Pro233. The length of its backbone skeleton also leads to forming two additional intramolecular bonds with the oxygen atom at C2 and the fluorine at C2'.

By examining the binding mode of Ro12-6377 at the / interface (Figs. 2d and 3d and Table 7), both the diazepine and phenyl (C) are positioned in front of the pocket. The methyl group attached at N1 forms two types of interactions. The first is a moderate H-donor bond which is given to the carbonyl of γ2Ser280. The seconds are three hydrophobic bond type Alkyl-Alkyl formed with side chains of γ2Ala283, γ2Arg284, and γ2Val300. Simultaneously, γ2Asp297 established two hydrogen bonds with the diazepine ring (A); strong H-donor has given to the oxygen attached at C2, and moderate H-acceptor bond received from C3. The phenyl (C) was involved in two hydrophobic interactions; type Pi-PiT-shaped with side chain of γ2Phe304, and type Pi-Alkyl with the side chain of β2Pro228. Its fluorine atom was oriented toward TM1:β2subunit, which resulted in halogen interaction with the lone pair of the oxygen situated in the main chain of Leu223, and in moderate H-acceptor bond forms with the side chain of Pro228. On the other hand, the phenyl (B) was inserted deeper in the binding pocket as its π-electron cloud shared Pi-Alkyl interaction with side chain of β2Pro228. Its NHCONHCH3 group attached at C7 was pointed towards the TM2:β2 helix. In addition to the two Alkyl-Alkyl interactions shared with side chains of Pro228 and Ile264 ofβ2 subunit, the NHCONHCH3 group also received one moderate H-bond from side chain of γ2Thr277 and simultaneously gave four moderate H-bonds, three of which to β2Gln224, and thus, leads to the formation of two intermolecular pentameric rings and one butameric ring contribute, as in the classical site, to the stability of Ro12-6377 at the binding site.

Table 7

Detailed binding interactions resulting from the molecular docking of co-crystallized DZP, Ro12-6377, and proflazepam (PLZ) at the TMD / interface

TMD γ2+/β2- interface
Ligand Name	Posita	Category	Type of interactions	From	From chemistry	To	To chemistry	Distb (A°)
DZP	3	Hydrogen bond	C-H bond	E:DZP403:H171	H-donor	E:SER280:OG	H-acceptor	2.77
	1			E:DZP403:H202		E:THR281:OG1		2.80
	Ring C	Hydrophobic	Pi-Pi T-shaped	E:PHE304	Pi-orbitals	E:DZP403	Pi-orbitals	5.81
	1		Alkyl	E:DZP403:C20	Alkyl	A:PRO228	Alkyl	4.25
	7			E:DZP403:CL		A:LEU223		4.64
				E:DZP403:CL		E:VAL300		4.38
	Ring B		Pi-Alkyl	E:DZP403	Pi-orbitals	E:VAL300		4.92
Ro12-6377	7	Hydrogen bond	Conventional H-bond	:*0:H	H-donor	A:GLN224:O	H-acceptor	2.43
	2’		C-H bond	A:PRO228:HD3		:*0:F		2.61
	7			E:THR277:HB		:*0:O		2.73
	2			E:ASP297:HA		:*0:O		2.28
	3			:*0:H		E:ASP297:OD1		2.89
	1			:*0:H		E:SER280:O		2.71
	7			:*0:H		A:GLN224:O		2.97
				:*0:H		E:THR281:OG1		2.78
				:*0:H		A:GLN224:O		2.94
	2’	Halogen	Halogen (Fluorine)	A:LEU223:O	Halogen acceptor	:*0:F	Halogen	2.56
	Ring C	Hydrophobic	Pi-Pi T-shaped	E:PHE304	Pi-orbitals	:*0	Pi-orbitals	5.40
	1		Alkyl	E:ALA283	Alkyl	:*0:C	Alkyl	4.28
				:*0:C		E:ARG284		3.96
				:*0:C		E:VAL300		3.90
	7			:*0:C		A:PRO228		4.75
				:*0:C		A:ILE264		3.81
	Ring B		Pi-Alkyl	:*0	Pi-orbitals	A:PRO228		5.11
	Ring C			:*0		A:PRO228		5.29
PLZ	1	Hydrogen bond	Conventional H-bond	A:GLN224:HE21	H-donor	:*0:O	H-acceptor	2.32
				A:GLN224:HE22		:*0:O		2.59
				:*0:H		E:THR281:OG1		1.92
	2’		C-H bond	A:PRO228:HD3		:*0:F		2.80
	2			E:THR277:HA		:*0:O		2.79
	7			E:ARG284:HD3		:*0:CL		2.83
	1			:*0:H		E:SER280:OG		2.74
				:*0:H		A:GLN224:O		2.81
				:*0:H		E:THR281:OG1		2.37
	7	Halogen	Halogen (Cl, Br, I)	E:ASP297:OD1	Halogen acceptor	:*0:CL	Halogen	3.18
	2’		Halogen (Fluorine)	A:LEU223:O		:*0:F		2.90
	Ring C	Hydrophobic	Pi-Pi T-shaped	E:PHE304	Pi-orbitals	:*0	Pi-orbitals	5.82
	7		Alkyl	:*0:CL	Alkyl	A:LEU223	Alkyl	4.29
	Ring B		Pi-Alkyl	:*0	Pi-orbitals	E:VAL300		4.92
	Ring C			:*0		A:LEU231		5.35
	1	Unfavorable	Unfavorable Acceptor-Acceptor	A:GLN224:O	H-acceptor	:*0:O	H-acceptor	2.83

*0: ligand structure

aPosition of interaction

bdistance of interaction

cinternal interaction

Dist

Proflazepam adopted a similar binding mode of DZP in the binding locus (Fig. S2 (d), Supplementary material). The binding mode of DZP at the / interface was previously discussed in detail by Kim et al. [7]. The uncommon binding interactions observed between the both modulators are related to the presence of two distinct binding groups on the structure of proflazepam: the fluorine atom at C2' and the long backbone skeleton of CH2CHOHCH2OH at N1. The phenyl (C) is placed in front of the pocket (Fig. 2d and 4d and Table 7), its π-electron cloud establishes two interactions of types Pi-PiT-shaped and Pi-Alkyl with side chains of γ2Phe304 and β2Leu231, respectively. Its fluorine atom points towards the TM1:β2subunit, where it receives a moderate H-bond from the side chain of Pro228 and establishes halogen interaction with the main oxygen of Leu223. The phenyl (B) binds at a higher level than rings (A) and (C), with its chlorine atom pointing towards the TM3:γ2 subunit. This orientation results in Pi-Alkyl interaction between the π-electron cloud and γ2Val300, and also leads the chlorine to accept the moderate H-bond from γ2Arg284, forms halogen interaction with the side chain of γ2Asp297, and interacts hydrophobically with β2Leu223. The binding modes of the oxygen atom attached at C2 in the three /,, and interfaces reflect the insensitivity of Ro12-6377 and proflazepam to participate in any interaction with the neighboring residues using this position. Otherwise, the / interface reflects the contribution of the oxygen atom at this position to enable the receptor potentiation by accepting, respectively, strong and moderate H-bonds from Asp297 and Thr277 of subunit γ2. The CH2CHOHCH2OH group penetrates deep into the binding site. Its pose and orientation towards the β2:A:TM2 helix are identical to that observed for the NHCONHCH3 group of Ro12-6377. The length of its backbone skeleton allowed it to simultaneously influence the γ2:TM2helix and the β2:A:TM1 by forming six H-bonds with γ2:TM2:Ser280, γ2:TM2:Thr281, and β2:A:TM1:Gln224.

The detailed mechanisms of interactions by which proflazepam modulates the two TMD and sites are appearing to be complementary to each other. As can be seen, the interface has been predicted to specifically modulate by the first hydroxyl group of CH2CHOHCH2OH, C3, the chlorine atom at C7, and the π-electron clouds of both phenyls (B) and (C), whereas the interface has been predicted to specifically modulate by the second hydroxyl group of CH2CHOHCH2OH, N4, the fluorine atom at C2', and the π-electron cloud of both phenyls (B) and (C). This selectivity in binding interactions between the two orthosteric binding sites allows us to hypothesize that they may be functional in a complementary manner, as previously observed for the two agonist (GABA) binding sites which we do not yet know an explanation for why they are structurally identical and functionally not equivalent [8]. Furthermore, the bulky structure of proflazepam allowed it to directly induce influence on the pore-lining helices β2:A:TM2, β2:C:TM2, and α1:D:TM2 by creating H-bonds with β2:A:Arg269,β2:A:Asn265, β2:C:Asn265, and interacting hydrophobically with α1:D:Leu269. Likewise, three hydrophobic interactions with β2:A:TM2:Met261, β2:C:TM2:Met261, and α1:D:TM2:Leu269 are observed for the DZP, whereas just one hydrophobic interaction with β2:A:TM2:Val258 is established for Ro12-6377. Evidently, rings (C) of both proflazepam and DZP adopt a similar hydrophobic interaction with α1:D:TM2:Leu269.

Obviously, at / interface, the binding mode of Ro12-6377 is the most influential on pore-lining by sharing five interactions with the residues of γ2:TM2helix and one interaction with β2:A:TM2 helix. Likewise, the binding mode of proflazepam is connected to the pore-lining by participating in five interactions with γ2:TM2 helix. Accordingly, both modulators exhibit common interactions with γ2:TM2residues: Thr281, Thr277, Arg284, and Ser280. Otherwise, by examining the binding interactions of DZP, its structure predicted to enrich the skeleton ofγ2:TM2 through two moderate H-bonds originating from the methyl group attached at N1 to Ser280 and from C3 to Thr281. These two interactions are identical between the three modulators.

The feature of interacting with the residues of the TM2 helices is of great importance as it leads both Ro12-6377 and proflazepam to directly induce motions in the chloride-channel lining, thereby, possibly contributing to the expansion of its diameter by opening the 9' gate by orienting the β2:C:Leu259 side chain towards one of the two adjacent α subunits. As mentioned earlier, this rotation is the main factor in the activation of the pLGICs family [7, 57].

Molecular dynamics simulation

MD simulation is used as a complementary tool to validate the docking results before they are approved in the drug-design process. MD simulation offers the peculiarity of treating biological systems as flexible entities. This flexibility allows for free integrations between the macromolecule binding site and the binding-ligand, resulting in binding modes (poses or interactions) that confirm or refute the results of molecular docking [58]. Or in some cases, it may lead to the release of the ligand from the binding site, and this is an undesirable defect especially if the ligand shows the highest stability in the docking simulation [59].

For these reasons, Ro12-6377 and proflazepam in complex with the four binding interfaces were subjected to MD simulations using the settled parameters cited in the “Materials and methods” section. Their dynamic behaviors were investigated by evaluating the response of the potential energy, U (Kcal/mol), over a time period of 1000 picoseconds (ps) (Fig. 5). During the first 500 ps of the simulation, the ECD / interface exhibits higher stability in complex with Ro12-6377 than in complex with proflazepam. Thereafter, from 500 ps until the end of the simulation, both complexes tend to have equivalent stability (6270.7 kcal/mol for Ro12-6377 and 6226.86 kcal/mol for proflazepam) (Fig. 5a). Moreover, within the three TMD binding interfaces, the two modulators exhibited stability equivalence throughout all the simulation periods (Fig. 5b–d).

Fig. 5

Evaluation the response of potential energy, U (kcal/mol), as function of time, t(ps), for Ro12-6377 and proflazepam in complex with a ECD /, b TMD /, c TMD /, and d TMD / interfaces

Later, deep analyses of binding modes, binding orientations, and binding interactions of Ro12-6377 and proflazepam within the four binding interfaces were performed and discussed between the two simulations.

MD simulation analysis of the classical binding site

At ECD / interface, the binding poses of Ro12-6377 are equivalents for the two simulations (Figs. 6a and 7a and Table 8). However, significant differences between the two binding orientations were detected, notably in substitutes at C2, C5, and C7. The new binding orientation predicted for phenyl (C) is driving the C2' bound fluorine atom to move away from the diazepine ring (A), which leads to the disappearance of the intramolecular interaction formed with C2. All the H-bonds formed with α1Ser206, γ2Phe77, and α1His102 were vanished and replaced by H-bonds given from side chain of α1Lys156 to the oxygen atom at C2, and from the NHCONHCH3 group at C7 to both α1Tyr160 and γ2Asp56. Similarly, the hydrophobic interactions suggested with Phe100 and Val203 of α1 subunit were also replaced by hydrophobic interactions created between the π-orbitals of Tyr160 and Tyr210 of the same subunit and the π-electron cloud of phenyl (C). In contrast, the hydrophobic interactions established with the subunit γ2 residues (Tyr58, Phe77, and Ala79) were preserved as the same.

Fig. 6

Binding modes resulting from the molecular dynamics simulation of Ro12-6377 at the interfaces of a ECD /, b TMD /, c TMD /, and d TMD / interfaces

Fig.7

Binding interactions resulting from the molecular dynamics simulation of Ro12-6377 atthe interfaces of a ECD /, b TMD /, c TMD /, and d TMD / interfaces

Table 8

Detailed binding interactions resulting from the molecular dynamics simulation Ro12-6377 and proflazepam (PLZ) at the classical site

ECD α1+/γ2- interface
Ligand Name	Posita	Category	Type of interactions	From	From chemistry	To	To chemistry	Distb (A°)
Ro12-6377	2	Hydrogen bond	Conventional hydrogen bond	D:LYS156:HZ3	H-donor	:*0:O	H-acceptor	1.67
	7			:*0:H		D:TYR160:OH		2.23
				:*0:H		D:TYR160:OH		2.27
			C-H bond	:*0:H		E:ASP56:OD1		2.63
	Ring C	Hydrophobic	Pi-Pi stacked	D:TYR160	Pi-orbitals	:*0	Pi-orbitals	4.46
	Ring B			E:TYR58		:*0		4.29
	Ring C		Pi-Pi T-shaped	D:TYR210		:*0		5.39
	Ring B			E:PHE77		:*0		4.68
	7		Alkyl	E:ALA79	Alkyl	:*0:C	Alkyl	3.59
	1		Pi-Alkyl	E:TYR58	Pi-orbitals	:*0:C	Alkyl	4.47
PLZ	4	Hydrogen bond	Conventional hydrogen bond	D:LYS156:HZ2	H-donor	:*0:N	H-acceptor	1.57
	1			:*0:H		D:SER159:O		1.79
	2’		C-H bond	D:LYS156:HE2		:*0:F		2.51
	Internalc			:*0:H		:*0:O		2.51
	Ring C	Electrostatic	Pi-Cation	D:LYS156:NZ	Positive	:*0	Pi-orbitals	4.25
	1	Hydrogen bond	Pi-Donor hydrogen bond	:*0:H	H-donor	D:TYR210		2.69
	Ring B	Hydrophobic	Pi-Alkyl	:*0	Pi-orbitals	D:VAL203	Alkyl	4.72

*0: ligand structure

aPosition of interaction

bdistance of interaction

cinternal interaction

The binding modes of proflazepam are inconsistent in the binding poses for the two simulations. The adequate binding pose generated by MD simulation was placed less deeply within the binding locus than that generated by molecular docking simulation (Fig. 8a). Unlike molecular docking simulation, this binding position emerged its structure to react as donor and acceptor with neighboring residues. The most notable differences in the binding orientations appear in the substitutes groups at N1,C5, and C7. By examining the binding interactions (Fig. 9a and Table 8), the moderate H-bond established between the first hydroxyl group of CH2CHOHCH2OH and the side chain of α1Tyr210 was preserved between the two simulations, whereas the remaining interactions were completely vanished and replaced by five interactions with the residues of the same subunit: side chain of Lys156 gives two H-donors to N4 and the fluorine atom bond at C2', and were also involved in Pi-Cation interaction with the π-electron cloud of phenyl (C). The main chain of Ser159 receives one H-acceptor bond from the second hydroxyl group of CH2CHOHCH2OH. A moderate intramolecular H-bond was formed between the oxygen atom at C2 and the CH2CHOHCH2OH group. Finally, hydrophobic interaction type Pi-Alkyl was observed between the π-electron cloud of phenyl (B) and the side chain of Val203.

Fig. 8

Binding modes resulting from the molecular dynamics simulation of proflazepam at the interfaces of a ECD /, b TMD /, c TMD /, and d TMD / interfaces

Fig. 9

Binding interactions resulting from the molecular dynamics simulation of proflazepam at the interfaces of a ECD /, b TMD /, c TMD /, and d TMD / interface

MD simulation analysis of the three TMD binding sites

The binding modes of Ro12-6377 for both TMD and interfaces are consistent in docking poses for the two simulations (Fig. 6b, c). However, notable differences were observed between the binding orientations of C5 and C7 substitutes. The binding interactions at both interfaces exhibit the tendency of Ro12-6377 to enhance its interactions with the residues of β2 subunit relative to those established with subunit α1(Figs. 7b, c and Tables 9 and 10). Briefly, at the interface, all interactions created with residues of subunit β2:A were preserved and stabilized by removing the moderate H-bond and the Pi-Alkyl interaction linking, respectively, the substitutes at N1 and C7 with the side chains of Leu285 and Phe289, and creating, instead of them, new moderate H-bond oriented from the methyl group at N1 towards the main carbonyl group of Asp282. In contrast, only interactions with two residues from α1:B subunit (Met236 and Thr237) were preserved. The Pi-Alkyl interaction bond of the side chain of Met236 to phenyl (C) was vanished and replaced with two H-bonds oriented from C7 substitute towards the main carbonyl group, and Pi-Sulfur interaction given from the sulfur atom in the side chain to phenyl (C). At the interface, the binding interactions with Thr237, Leu240, and Leu232 of α1:D subunit were vanished and replaced by new interactions with residues from the β2:C subunit: Arg269, Leu285, and Val258. Side chains of both Leu285 and Val258 were, respectively, involved in Alkyl-Alkyl interactions with the substitutes at N1 and C7, while the side chain of Arg269 established a strong H-bond with the oxygen atom at C2. Likewise, the intramolecular interaction observed between N1 and the fluorine atom at C2' also vanished and created, instead of it, additional interaction between this latter and the main carbonyl of Met286. At the α1:D subunit, likewise to B:Met236, D:Met236 was replacing the Pi-Alkyl interaction connected its side chain to phenyl (C) by two moderate H-bonds and one Pi-Sulfur interaction. The two moderate H-bonds are oriented from the NHCONHCH3 group at C7 towards the main carbonyl group, whereas the Pi-Sulfur interaction appears between the sulfur atom in the side chain and the phenyl (B).

Table 9

Detailed binding interactions resulting from the molecular dynamics simulation of Ro12-6377 and proflazepam (PLZ) at the TMD / interface

TMD β2+(A)/α1-(B) interface
Ligand Name	Posita	Category	Type of interactions	From	From chemistry	To	To chemistry	Distb (A°)
Ro12-6377	7	Hydrogen bond	Conventional hydrogen bond	:*0:H	H-donor	B:MET236:O	H-acceptor	2.51
				:*0:H		B:MET236:O		1.85
	1		C-H bond	:*0:H		A:ASP282:O		2.79
	7			:*0:H		B:THR237:OG1		3.06
	Ring C	Other	Pi-Sulfur	B:MET236:SD	Sulfur	:*0	Pi-orbitals	3.90
	Ring B	Hydrophobic	Pi-Pi stacked	A:PHE289	Pi-orbitals	:*0		4.62
	1		Alkyl	:*0:C	Alkyl	A:LEU285	Alkyl	5.05
				:*0:C		A:MET286		4.81
	7			:*0:C		A:VAL258		5.23
	Ring C		Pi-Alkyl	:*0	Pi-orbitals	A:MET286		4.84
PLZ	Internalc	Hydrogen bond	Conventional hydrogen bond	:*0:H	H-donor	:*0:O	H-acceptor	1.63
	1			:*0:H		B:ILE228:O		1.85
			C-H bond	A:ARG269:HD2		:*0:O		2.82
				:*0:H		A:ASN265:OD1		2.71
	Ring B	Hydrophobic	Pi-Pi stacked	A:PHE289	Pi-orbitals	:*0	Pi-orbitals	4.73
	7		Alkyl	:*0:CL	Alkyl	A:MET286	Alkyl	5.41
				:*0:CL		B:LEU232		4.88
				:*0:CL		B:PRO233		4.64
				:*0:CL		B:MET236		4.14
	7		Pi-Alkyl	A:PHE289	Pi-orbitals	:*0:CL		4.39
	Ring B			:*0		B:PRO233		4.39

*0: ligand structure

aPosition of interaction

bdistance of interaction

cinternal interaction

Table 10

Detailed binding interactions resulting from the molecular dynamics simulation of Ro12-6377 and proflazepam (PLZ) at the TMD / interface

TMD β2+(C)/α1-(D) interface
Ligand Name	Posita	Category	Type of interactions	From	From chemistry	To	To chemistry	Distb (A°)
Ro12-6377	2	Hydrogen bond	Conventional hydrogen bond	C:ARG269:HH12	H-donor	:*0:O	H-acceptor	1.62
	7			:*0:H	H-donor	D:MET236:O		2.55
	2’	Hydrogen bond; Halogen	C-H bond; halogen (Fluorine)	C:MET286:HA	H-donor; halogen acceptor	:*0:F	H-acceptor; halogen	2.53
	3	Hydrogen bond	C-H bond	:*0:H	H-donor	D:ILE228:O	H-acceptor	3.03
				:*0:H		D:ILE228:O		2.75
	7			:*0:H		D:MET236:O		2.81
	2’	Halogen	Halogen (Fluorine)	C:MET286:O	Halogen acceptor	:*0:F	Halogen	3.32
	Ring B	Other	Pi-Sulfur	D:MET236:SD	Sulfur	:*0	Pi-orbitals	5.26
		Hydrophobic	Pi-Pi stacked	C:PHE289	Pi-orbitals	:*0		4.95
	1		Alkyl	:*0:C	Alkyl	C:LEU285	Alkyl	4.68
				:*0:C		C:MET286		4.54
	7			:*0:C		C:VAL258		4.89
	Ring B		Pi-Alkyl	:*0	Pi-orbitals	D:PRO233		4.88
	Ring C			:*0		C:MET286		5.12
PLZ	2	Hydrogen bond	Conventional hydrogen bond	C:ASN265:HD22	H-donor	:*0:O	H-acceptor	2.20
	1			C:ARG269:HH12		:*0:O		1.42
	Internalc	Hydrogen bond; halogen	Conventional hydrogen bond; halogen (Fluorine)	:*0:H	H-donor; halogen acceptor	:*0:F	H-acceptor; halogen	1.81
	1	Hydrogen bond	Conventional hydrogen bond	:*0:H	H-donor	D:ILE228:O	H-acceptor	1.65
	Internalc		C-H bond	:*0:H	H-donor	:*0:F		2.50
	Internalc	Halogen	Halogen (Fluorine)	:*0:N	Halogen acceptor	:*0:F	Halogen	3.44
	Ring B	Other	Pi-Sulfur	D:MET236:SD	Sulfur	:*0	Pi-orbitals	5.49
		Hydrophobic	Pi-Pi stacked	C:PHE289	Pi-orbitals	:*0		4.64
	7		Alkyl	:*0:CL	Alkyl	D:PRO233	Alkyl	4.81
				:*0:CL		D:MET236		5.26
			Pi-Alkyl	C:PHE289	Pi-orbitals	:*0:CL		4.94
	Ring B			:*0		D:PRO233		4.44
	Ring C			:*0		D:PRO233		5.31
				:*0		D:LEU269		5.36

*0: ligand structure

aPosition of interaction

bdistance of interaction

cinternal interaction

The binding modes of proflazepam for the TMD interface are inconsistent in docking poses for the two simulations (Fig. 8b). Accordingly, the differences between the two binding orientations were detected for the entire structure of proflazepam. The adequate binding pose resulting from the MD simulation was inserted more deeply into the binding locus in such a way that the fused benzodiazepine rings and the substitute at N1 have facial alignment with the β2:TM3 helix. The binding interactions after the MD simulation (Fig. 9b and Table 9) suggested the lack of the two H-bonds established with β2Leu285 and β2Arg269, as well as the three Pi-Alkyl interactions which were created with β2Met286, α1Leu232, and α1Pro233. Alternatively, α1Ile228, β2Phe289, α1Pro233, and β2Met286 were estimated to participate in six interactions with phenyl (B) and the substitutes at N1 and C7. The second hydroxyl group of CH2CHOHCH2OH gave a strong H-bond to the main carboxyl function of α1Ile228. The chlorine at C7 was involved in three hydrophobic interactions, two are of the Alkyl-Alkyl type created with side chains of β2Met286 and α1Pro233 and one is of the Pi-Alkyl type created with the Pi-orbitals of β2Phe289. Finally, the Pi-orbitals of phenyl (B) have participated through Pi-Pi stacked and Pi-Alkyl interactions with the side chains of β2Phe289 and α1Pro233, respectively.

The binding modes of proflazepam for the TMD interface are equivalents in the docking poses for the two simulations (Fig. 8c). The most pronounced differences in the binding orientations have appeared for the substitutes groups at N1, C5, and C7. The molecular docking and MD simulations overlap in that the residues β2Asn265, β2Phe289, α1Ile228, α1Pro233, and α1Leu269 are essential parts of the modulation by proflazepam (Fig. 9c and Table 10). However, this does not necessarily mean that they have retained the same nature and orientation of their binding interactions shared with the proflazepam structure for both simulations, which is an observation common to all the binding sites studied in this paper. The binding pose after the MD simulation was stabilized by three intramolecular interactions orienting from the CH2CHOHCH2OH group and N1 towards the fluorine atom bond at C2'. The four hydrophobic interactions created for both phenyls (B) and (C) and the H-bond linked α1Ile 228 to the CH2CHOHCH2OH group were preserved as the same. β2Asn265 reduced the two H-bonds established with N4 and the fluorine atom at C2' to a single strong H-bond directed from its amide group towards the bound oxygen atom at C2. α1Ile228 loses the halogen interaction established with the fluorine atom at C2' and conserves the H-bond accepted from the CH2CHOHCH2OH group. α1Pro233 tends to be involved in hydrophobic interactions using its side chain. Thus, in addition to the two Pi-Alkyl interactions created with phenyl (B) and phenyl (C), it prefers to create an Alkyl-Alkyl interaction with the chlorine atom at C7 rather than the strong H-bond created with the substitute at N1. β2Phe289 shows an additional hydrophobic interaction type Pi-Alkyl with the chlorine atom at C7. Finally, three interactions with two new residues were observed for phenyl (B) and the substitutes at N1 and C7. The second hydroxyl group of CH2CHOHCH2OH accepts a strong H-bond from the side chain of β2Arg269, and both the π-electron cloud of phenyl (B) and the chlorine at C7 formed, respectively, Pi-Sulfur and Alkyl-Alkyl interactions with the side chain of α1Met236.

At TMD / interface, the binding modes of Ro12-6377 show significant differences for the two simulations. The adequate binding mode generated by MD simulation was inserted more deeply into the binding locus, so that phenyl (B) is positioned between β2:TM1 and γ2:TM2 helices and phenyl (C) has facial alignment with the β2:TM2 and γ2:TM3 helices (Fig. 6d). MD simulation reduced the number of interactions suggested by molecular docking to less than half (Fig. 7d and Table 11), or else, suggested three interactions with γ2Lys285 and β2Leu268 instead of those established with β2Pro228, β2Leu223, γ2Thr277, γ2Thr281, γ2Phe304, and γ2Ala283. Both γ2Lys285 and β2Leu268 interact with the substitution group at C7 through two H-bonds and one Alkyl-Alkyl interaction orienting, respectively, from their side chains towards the lone pairs of the oxygen atom and the methyl group. On the other hand, the number of interactions created with γ2Asp297 and β2Gln224 is reduced by the factor of one H-bond for each. In contrast, interactions with γ2Ser280, γ2Arg284, γ2Val300, and β2Ile264 residues are kept identical between the two simulations.

Table 11

Detailed binding interactions resulting from the molecular dynamics simulation of Ro12-6377 and proflazepam (PLZ) at the TMD / interface

TMD γ2+/β2- interface
Ligand Name	Posita	Category	Type of interactions	From	From chemistry	To	To chemistry	Distb (A°)
Ro12-6377	7	Hydrogen bond	Conventional hydrogen bond	E:LYS285:HZ2	H-donor	:*0:O	H-acceptor	2.09
				:*0:H		A:GLN224:O		2.94
				:*0:H		A:GLN224:O		1.74
			C-H bond	E:LYS285:HE2		:*0:O		2.73
	2			E:ASP297:HA		:*0:O		2.40
	1			:*0:H		E:SER280:O		2.96
		Hydrophobic	Alkyl	:*0:C	Alkyl	E:ARG284	Alkyl	4.71
				:*0:C		E:VAL300		4.13
	7			:*0:C		A:ILE264		4.75
				:*0:C		A:LEU268		4.66
PLZ	1	Hydrogen bond	Conventional hydrogen bond	:*0:H	H-donor	E:THR281:OG1	H-acceptor	1.81
				:*0:H		A:GLN224:O		2.06
			C-H bond	:*0:H		E:THR277:O		2.86
				:*0:H		E:THR277:O		2.94
				:*0:H		E:SER280:OG		2.51
				:*0:H		E:THR281:OG1		2.60
	Ring C	Other	Pi-Sulfur	A:MET227:SD	Sulfur	:*0	Pi-orbitals	5.18
	7	Hydrophobic	Alkyl	:*0:CL	Alkyl	A:LEU223	Alkyl	4.37
				:*0:CL		E:VAL300		4.07
	Ring B		Pi-Alkyl	:*0	Pi-orbitals	E:VAL300		5.15

*0: ligand structure

aPosition of interaction

bdistance of interaction

Unlike Ro12-6377, the binding modes of proflazepam for the TMD  interface are equivalents for the two simulations (Fig. 8d). However, notable differences in the binding orientations of the substitution groups at C5 and C7 have appeared. As previously noted for Ro12-6377, after MD simulation, the number of interactions was almost reduced to half (Fig. 9d and Table 11). All the interactions established with β2Pro228, γ2Arg284, γ2Asp297, β2Leu231, and γ2Phe304 have vanished, and instead of them, Pi-Sulfur interaction was observed between the sulfur atom at β2Met227 side chain and the Pi-orbitals of phenyl (C). The three H-bonds linked β2Gln224 to the CH2CHOHCH2OH group were reduced into one moderate H-acceptor bond oriented from the second hydroxyl group towards the main carbonyl group of β2Gln224. The halogen interaction established between β2Leu223 and the fluorine at C2' vanished, and the H-bond connected C2 to γ2Thr277 was replaced by two moderate H-bonds given from the CH2CHOHCH2OH group to the main carbonyl group of the same residue. γ2Val300 was involved in additional hydrophobic interaction type Alkyl-Alkyl with the chlorine atom at C7. Finally, interactions with γ2Thr281 and γ2Ser280 residues are kept identical between the two simulations.

Overall, the MD results do not agree with those observed in the molecular docking analysis that showed interactions with α1His102. As mentioned earlier, this residue ensures the recognition of classical BDZ by the ECD interface. The MD results coincided well with those noted in the molecular docking analysis that indicated the features of Ro12-6377 and proflazepam to directly connect the pore-lining residues. The bulky structure of Ro12-6377 and proflazepam is the key factor in the deep penetration towards the TM2 helices, in particular, the long backbone skeleton of the NHCONHCH3 and CH2CHOHCH2OH groups, where most interactions with the TM2 helices have been observed. At TMD , the MD simulation generated two interactions for Ro12-6377 with β2:C:TM2:Arg269 and β2:C:TM2:Val 258 residues. In addition, at TMD / interface, three interactions were generated with β2:A:TM2:Leu268 and γ2:TM2:Lys285 instead of those previously observed with TM2:Thr277, TM2:Thr281, and TM2:Ala283of subunitγ2. Likewise, at TMD interface, it suggested an additional interaction for the proflazepam with β2:C:TM2:Arg269, and at TMD  interface, suggested a lack of interactions with γ2:TM2:Arg 284.

Pharmacokinetic and drug-likeness prediction

As shown in Table S2, Supplementary materials, all BDZ compounds are estimated to have high gastrointestinal absorption and nearly all are able to cross the blood–brain barrier. Almost all BDZ molecules are not affected by p-gp efflux pump. Besides, all tested compounds respect Lipinski, Veber, Egan, Ghose, and Muegge rules.

Moreover, Ro12-6377 and proflazepam can successfully penetrate the blood–brain barrier. They are also estimated to be active effluxes by P-glycoprotein transporter and to act as non-inhibitors towards CYPisoenzymes, except for CYP2D6 in proflazepam. Interestingly, these two compounds share one favorable characteristic in which they do not inhibit CYP2C19 and CYP3A4 enzymes that might be responsible for the hepatic clearance or the formation of active metabolites of BDZ. Indeed, it is well established that the BDZs are primarily metabolized via CYP2C19 and CYP3A4 and the inhibition of them can cause drug-drug interactions (DDIs) [60, 61].

The QSAR study was carried out on the BDZ data set previously investigated through molecular docking simulation (Table 1). The data set includes a total of 65 compounds that were used to generate the PLS regression model and evaluate its performance. Using fivefold cross-validation, we randomly split the 65 observations into two sets, a training set containing 52 of the data points, and a test set containing the remaining 13 observations. The molecular descriptors are coded into the term of variables: x1 = qN1, x2 = qC3, x3 = qN4, x4 = qC6, x5 = qC7, x6 = qC8, x7 = qC9, x8 = qC2', x9 = qC6', x10 = πC7, x11 = πC2', x12 = HA, x13 = HD, x14 = DM, x15 = Log P, x16 = MR, x17 = FT. The numerical values are summarized in Table S3, supplementary materials.

Table 12 reports the observed (ȳ) and predicted (ŷ) biological activities with the corresponding studentized deleted residual and the leverage (h) values of the training and test set compounds of a sample among the 10,000 simulations generated in this study. In this case, six outliers (marked in bold) were detected for the activity response. Four observations have greater than the threshold |2| and three observations with large h values. Consequently, the training and test sets were reduced to 47 and 12 observations, respectively.

Table 12

Studentized deleted residual values () and the leverage values (h) of the training and test set compounds

CompN°	Yinorm	Yhat	ri∗	hii	CompN°	Yinorm	Yhat	ri∗	hii
Training set
1	8.2954	8.1891	0.1713	0.2220	35	8.2954	8.5488	− 0.4167	0.2493
3	8.4000	9.1331	− 1.1372	0.1352	36	9.2113	9.1223	0.1354	0.1266
4	8.4393	9.7341	− 2.0056	0.0804	37	9.1982	9.5757	− 0.5616	0.0800
6	9.3159	9.8727	− 0.9261	0.2553	39	8.3346	8.8803	− 0.8621	0.1764
7	9.6430	9.2510	0.5873	0.0922	40	8.1253	8.9133	− 1.2514	0.1698
8	9.6823	10.0519	− 0.5490	0.0775	43	8.9627	9.6633	− 1.4642	0.5146
9	9.7346	9.7370	− 0.0036	0.1082	44	8.5309	9.4693	− 1.5789	0.2455
10	9.8393	9.5026	0.4970	0.0669	45	11.0692	10.9866	0.1345	0.2361
11	9.9701	9.8377	0.1953	0.0698	47	10.6636	9.8739	1.2419	0.1540
13	10.1271	9.7364	0.5867	0.0962	48	9.5645	10.1277	− 0.8436	0.0843
14	10.3234	10.4628	− 0.2098	0.1072	49	10.1926	9.7383	0.6890	0.1121
15	10.5851	9.9585	0.9361	0.0763	50	10.8206	8.6790	4.5068	0.3328
16	10.8075	11.0797	− 0.4154	0.1287	51	8.8580	9.5543	− 1.0520	0.0922
17	10.8206	10.0065	1.2313	0.0859	52	10.1010	10.2745	− 0.2602	0.1000
19	11.0561	10.9483	0.1691	0.1788	53	10.2711	9.9006	0.5519	0.0827
20	11.3832	10.6639	1.0698	0.0625	54	11.0169	10.7783	0.3697	0.1553
21	11.4356	11.6426	− 0.3182	0.1425	55	8.4916	9.3474	− 1.2840	0.0683
22	11.4356	11.1475	0.4431	0.1420	56	8.9234	9.4633	− 0.7967	0.0581
23	11.5402	10.0976	2.2560	0.0778	58	10.1010	9.8097	0.4441	0.1263
25	10.1664	10.7639	− 1.0675	0.3503	59	9.3421	9.2794	0.0929	0.0785
27	10.5066	10.2804	0.3321	0.0600	60	8.5309	9.3314	− 1.2208	0.1014
29	10.9776	10.8453	0.2004	0.1183	61	9.6823	9.8128	− 0.1930	0.0745
30	11.1477	11.5546	− 0.6271	0.1413	62	11.0430	10.4101	0.9384	0.0621
31	10.6636	11.4656	− 1.2521	0.1410	63	10.2711	10.0263	0.3632	0.0791
32	10.4542	9.7895	1.0285	0.1353	64	9.7215	9.7755	− 0.0789	0.0547
33	11.3309	9.9877	2.2297	0.1836	65	9.5253	9.8470	− 0.4706	0.0503
Test set
2	8.5167	8.6018	− 1.2800	0.3654	34	9.5531	9.1409	0.3076	0.3440
5	9.3272	8.8671	0.5421	0.3889	38	9.2475	9.0331	− 0.5152	0.9222
12	10.2041	10.2062	− 0.3840	0.5716	41	10.1111	10.3151	− 2.3002	0.4636
18	11.1873	9.7623	1.8961	0.5383	42	9.7922	9.4752	0.8742	0.8557
24	11.0146	10.1997	0.7054	0.3238	46	11.8516	12.0887	− 1.4227	0.9560
26	10.8020	10.4640	− 0.1934	0.2681	57	10.0314	10.2422	− 0.7454	0.4415
28	11.6789	10.8360	1.6756	0.5610

Figures 10 and 11 show the results of the 10,000 simulations obtained for the selected biological activity. According to Fig. 10, the best subsets of variables are those with the highest probability of occurrence according to the BIC criterion. As a result, the best subset of variables to select is that containing six variables. The intercept is systematically included whatever the model. The next step is determining the variables retained in the subset. According to Fig. 11, it appears that the six best variables are x3, x11, x14, x10, x15, and x16. This ranking order is a function of their probability of occurrence. By combining the results of Figs. 10 and 11, the best variable subset contains: qN4, πC2', DM, πC7, log P, and MR (Table 13).

Fig. 10

Box plots of the distribution of the Bayesian Information Criterion (BIC) by number of molecular descriptors

Fig. 11

Probability of occurrence of selected molecular descriptors

Table 13

The optimal variables to generate the PLS model

CompN°	qN4	π7	π2'	MD (debye)	Log P	MR	CompN°	qN4	π7	π2'	MD (debye)	Log P	MR
1	− 0.608	− 1.230	0.000	6.301	1.643	79.708	34	− 0.665	− 1.230	0.140	6.893	1.786	79.924
2	− 0.618	− 1.230	0.000	6.583	1.779	76.592	35	− 0.655	− 1.030	0.140	4.302	1.743	93.660
3	− 0.590	− 0.570	0.000	2.121	2.328	80.729	36	− 0.573	0.710	0.000	4.083	4.031	85.262
4	− 0.620	0.000	0.000	5.061	2.608	71.891	37	− 0.607	0.710	0.000	3.104	3.304	87.392
5	− 0.655	− 1.340	0.140	6.721	2.123	82.211	38	− 0.584	0.710	0.000	3.049	3.857	91.754
6	− 0.698	− 1.230	0.710	6.791	2.383	81.396	39	− 0.569	− 0.280	0.000	1.785	2.475	87.181
7	− 0.620	− 0.650	0.000	3.137	2.320	78.475	40	− 0.587	0.710	0.000	3.222	4.130	93.618
8	− 0.613	0.140	0.000	3.369	2.751	72.108	41	− 0.642	0.710	0.140	3.078	2.529	86.321
9	− 0.619	1.020	0.000	5.662	3.566	81.533	42	− 0.640	0.710	0.140	5.437	1.843	98.979
10	− 0.634	− 0.570	0.140	3.218	2.471	80.945	43	− 0.634	0.710	0.140	5.140	1.898	92.283
11	− 0.631	0.820	0.000	4.936	3.345	81.578	44	− 0.657	− 0.280	0.710	4.125	4.070	99.938
12	− 0.657	0.000	0.140	5.752	2.751	72.108	45	− 0.736	0.710	0.140	3.539	3.923	81.406
13	− 0.670	− 0.550	0.140	3.274	2.308	82.510	46	− 0.766	− 0.280	0.710	2.404	3.721	87.510
14	− 0.614	0.880	0.000	1.909	3.486	77.865	47	− 0.698	− 0.280	0.140	3.063	3.123	86.038
15	− 0.601	0.710	0.000	3.228	3.076	79.812	48	− 0.660	0.710	0.000	2.958	3.645	84.306
16	− 0.679	0.710	0.710	3.462	4.284	89.421	49	− 0.609	0.710	0.000	4.751	2.787	81.010
17	− 0.659	0.460	0.140	4.057	4.085	82.395	50	− 0.541	0.710	0.140	6.653	4.124	84.764
18	− 0.622	− 0.280	0.140	3.209	2.555	81.544	51	− 0.647	0.560	0.000	5.677	3.635	81.974
19	− 0.673	− 0.280	0.880	3.692	3.426	84.185	52	− 0.674	0.000	0.140	5.429	2.893	72.324
20	− 0.665	0.710	0.140	4.140	3.355	76.912	53	− 0.648	0.000	0.140	5.768	2.615	75.224
21	− 0.696	0.710	0.710	3.938	3.816	81.501	54	− 0.666	0.000	0.710	5.557	3.076	79.812
22	− 0.696	− 0.280	0.710	2.488	3.152	83.016	55	− 0.605	0.000	0.000	4.500	3.212	76.696
23	− 0.649	− 0.280	0.140	2.378	2.691	78.428	56	− 0.635	0.000	0.140	5.583	3.219	80.028
24	− 0.637	0.140	0.140	4.373	2.757	75.440	57	− 0.684	0.000	0.140	4.474	3.497	77.129
25	− 0.597	0.140	0.000	3.338	2.615	57.224	58	− 0.676	0.000	0.140	6.296	3.264	77.149
26	− 0.654	0.140	0.140	4.242	2.893	72.324	59	− 0.620	0.000	0.140	5.954	3.219	80.028
27	− 0.624	0.710	0.000	3.257	3.212	76.696	60	− 0.632	0.000	0.140	4.116	3.823	84.833
28	− 0.689	0.710	0.140	3.645	3.497	77.129	61	− 0.611	0.710	0.000	2.420	3.680	84.617
29	− 0.700	0.710	0.140	3.646	3.361	80.245	62	− 0.662	0.710	0.140	3.363	3.959	81.717
30	− 0.681	0.710	0.710	3.557	3.959	81.717	63	− 0.668	0.560	0.140	5.148	3.868	81.954
31	− 0.710	0.710	0.710	3.486	4.420	86.306	64	− 0.614	0.710	0.000	3.589	3.816	81.501
32	− 0.615	− 0.280	0.000	1.381	2.548	78.212	65	− 0.617	0.710	0.000	3.336	3.725	81.737
33	− 0.663	− 0.280	0.710	3.365	3.016	97.151

In order to highlight the weight of each molecular descriptor, the regression models are written with scaled variables. The standardized regression coefficient value of each descriptor highlights the relative importance of the descriptors in determination of biological activity of the compounds.

The final QSAR model with 95% confidence interval of the regression coefficient is:

is the biological activity of [3H] diazepam derivatives, where and . S: the standard deviation corresponding to the biological response. S: the standard deviation corresponding to the jth descriptor.

The best subset of molecular descriptors providing a good prediction of the response variable corresponds to: x3 = qN4, x10 = πC7, x11 = πC2', x14 = DM, x15 = Log P, x16 = MR.

Table 14 summarizes the statistical indicators used for internal and external validation. According to the goodness of fit statistics, 63.2% of the variability in BDZ activity around its mean is explained by the PLS regression equation. The quality of models can be judged and compared based on the values. The F-statistics reveal the significance of the PLS regression equations. The obtained p-value shows that the model is statistically highly significant at 95%. Moreover, it is well known that cross-validation is useful for overcoming the problem of overfitting [32]. This problem refers to a situation when the model requires more information than the data can provide. Indeed, in our case the difference between R2 and is much less than the threshold of 0.30, confirming that the PLS regression models are not overfitted. Additionally, our model exhibits a high value of (0.813). This result confirms that the resulting QSAR model has good external predictability and robustness.

Table 14

Quality and validation metrics

Model	Goodness of fit				Goodness of prediction
	R2	Radj2	F	p-value	Qloo2	QF32
Y (ntr = 52; nts = 13)	0.632	0.584	12.806	6.2050e − 07	0.639	0.813

ntr is the training set used for building the PLS regression equation and nts is the test set used to verify a model’s predictive ability for new untested molecules

Molecular descriptors with positive regression coefficients indicate positive correlations with observed activity. Therefore, their increase will improve the GABAA/BDZ response. Our model demonstrates a positive contribution from the hydrophobic constants of the substitutes at C7 and C2', as well as demonstrates a negative contribution from molecular lipophilicity (log P), molecular polarity (DM), molecular size (MR), and net charge of N4(qN4).

According to the standardized regression coefficient values, the hydrophobicity in positions 7 and 2’ are the main properties for determining the biological activities of the studied compounds. This result is in good agreement with previous reports that the 7 and 2’ positions are, respectively, the first and second key positions in the BDZ structures to influence binding with GABAARs [14]. Otherwise, the lipophilic behavior shows the least effect on activity. Molecular lipophilicity is directly influenced by both molecular hydrophobicity and polarity [62]. In our case, the positive effects of hydrophobicity at positions 7 and 2’ on the observed activity seem to be two times more significant than the negative effect of molecular polarity. The lipophilic feature plays a pivotal role in understanding the pharmacokinetics parameters, pharmacodynamics, and toxicological profile of drugs. Besides, it exhibits an important influence on host–guest interaction and drug binding affinity [46]. BDZ derivatives are relatively weak bases with highly lipophilic characters. At physiological pH, this lipophilicity explains their strong binding to plasma proteins (70–99%) as well as their high penetration through the blood–brain barrier [63].

N4 net charge and molar refractivity exhibit about twice the negative influence of molecular lipophilicity on activity. Reducing the negative charge of N4 is important to avoid the formation of water-unstable BDZ salts. The imine group (N4) is the most basic nitrogen in the classical BDZ structure. Since the amide group at positions 1 and 2 has a non-basic character, the lone pair of N4 can easily protonate when placed in a strongly acidic environment. Thus, it leads to the formation of BDZ salts (iminium ion). Unfortunately, the salts of strong acids are unstable in an aqueous medium and are therefore undergoing sequential hydrolysis of the imine and amide groups, respectively. The imine hydrolysis reaction is reversible. In contrast, hydrolysis of amide leads to the formation of inactive products and, consequently, eliminates activity towards GABAA receptors [64]. Most of the classical BDZ agents we study here provide unstable salts and are relatively water-insoluble than drugs formed from heterocyclic BDZs. In the latter, the amide group is protected in the form of heterocyclic groups such as ImidazoBDZs and TriazoloBDZs. Hence, hydrolysis of the amide does not occur and the reaction does not lead to inactive products. To improve the solubility of classicalBDZ salts in water-soluble injections, in addition to water, it is necessary to use co-solvents such as PEG 400, propylene glycol, 10% ethyl alcohol, and 2% benzyl alcohol [65].

Our results support those of Maddalena and So [14, 15] by asserting that increased hydrophobicity at position 7 is necessary to obtain highly potent BDZ analogs. Besides, it provides complementary insights on how hydrophobicity at C2’ position, the net charge of N4, and molar refractivity, polarity, and lipophilicity of the entire molecules would be ameliorated to achieve the optimal activity. In particular, the challenge is to provide structures meeting simultaneously the requirements of low lipophilicity/low polarity, since they are naturally anti-correlated. Here, our model suggests that the influence of molecular polarity on activity outweighs that of molecular lipophilicity by 7%. Accordingly, efforts should be devoted primarily to reducing molecular polarity rather than molecular lipophilicity.

After excluding outliers, a comprehensive analysis of the detailed binding interactions with the four binding interfaces was performed for the training and test data sets. Subsequently, the hydrophobic interactions with C7 and C2' substitutes and the electrostatic interactions with N4 were selected and collected in Tables S4 to S7, in supplementary materials. As can be seen, the richness of the C7 position by hydrophobic interactions is attributed in most cases to the presence of chlorine atoms that exhibited Alkyl-Alkyl or Pi-Alkyl interactions with neighboring residues. Significantly, the response of chlorine to establish hydrophobic interactions with the four binding interfaces is different between ECD and TMD. At ECD / interface, it tends to act as an acceptor of interactions, while at the three TMD interfaces it acts as a generator of interactions. C2' position shows a total lack of hydrophobic interactions with the four binding interfaces, except for the case of the CF3 group in Ro05-3590 and the chlorine atom in Ro05-4608. This deficiency maybe because most of the compounds of our data set contains small substituents (H and F) which mainly tend to interact electrostatically rather than hydrophobically. The N4 atom tends to react as a hydrogen-acceptor from neighboring residues. Importantly, in the ECD / interface, all interactions observed for N4 are received from α1His 102 side chain. As reported previously in the molecular docking section, this residue is known to be important in the recognition of classical BDZ. This result does not correlate with previous analyses that support cationic interactions at this position [66].

Conclusion

In this investigation, a combination of in silico approaches including molecular docking/dynamic simulations, and QSAR analysis have been performed to achieve two purposes: elucidate the binding mechanism by which a dataset of [3H]diazepam derivatives allosterically modulates GABAA receptor α1β2γ2 subtypes and identify the structural details that contribute to ameliorate the α1β2γ2/BDZ response. Examination of binding affinities revealed that the known ECD is the target for the majority of classical benzodiazepines. However, the tendency of the remainder to bind mainly on the binding interfaces included in TMD, or in some cases, to act on both ECD and TMD binding sites simultaneously, cannot be overlooked. This result opens the way for further studies that may combine binding in these binding sites with diversity in the activities of benzodiazepines. Binding affinity-based screening identified Ro12-6377and proflazepam as the best modulators for the four binding interfaces. By monitoring the dynamic behaviors over a time period of 1000 ps, the two modulators were observed to have equivalent stability within the four binding sites under study. Binding modes after MD simulation were altered from the structures generated by molecular docking. Thus, several differences in the binding interactions with key residues were detected between the two simulations. Importantly, interactions with pore-lining residues have been suggested for both modulators. The combination of ADMEprediction/Drug-likeness prediction shows their good pharmacokinetic properties as well as their compliance with all drug-likeness rules. Furthermore, the developed QSAR model yielded satisfactory statistical results that explain 63.2% of the variability in benzodiazepine activity. Its stability and predictive power were ensured based on internal and external validation indicators: =0.584, F = 12.806; p-value = 6.2050e − 07, =0.639, and =0.813. The model equation demonstrates a positive contribution from the hydrophobicity of the substitutes at C7 and C2', as well as demonstrates a negative contribution from molecular lipophilicity, molecular polarity, molecular size, and net charge of N4. The model results agree well with previous findings indicating that the increase in hydrophobicity at 7-position mainly contributes to the enhancement of BDZ activity. Finally, the combination of the results of both methods: QSAR and molecular docking, shows that the hydrophobic interactions at the 7-position are mostly attributed to the substitutions of chlorine atoms. The latter tends to act as an acceptor of interactions in the ECD binding interfaces and as a generator of interactions in the three TMD binding interfaces.

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 5163 KB)