{'fraction_non_alphanumeric': 0.037057657617531124, 'fraction_numerical': 0.026924704362737455, 'mean_word_length': 5.032240730789898, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Orthosteric and benzodiazepine cavities of the α1β2γ2 GABAA receptor: insights from experimentally validated in silico methods

γ-aminobutyric acid-type A (GABAA) receptors mediate fast synaptic inhibition in the central nervous system of mammals. They are modulated via several sites by numerous compounds, which include GABA, benzodiazepines, ethanol, neurosteroids and anaesthetics among others. Due to their potential as targets of novel drugs, a detailed knowledge of their structure–function relationships is needed. Here, we present the model of the α1β2γ2 subtype GABAA receptor in the APO state and in complex with selected ligands, including agonists, antagonists and allosteric modulators. The model is based on the crystallographic structure of the human β3 homopentamer GABAA receptor. The complexes were refined using atomistic molecular dynamics simulations. This allowed a broad description of the binding modes and the detection of important interactions in agreement with experimental information. From the best of our knowledge, this is the only model of the α1β2γ2 GABAA receptor that represents altogether the desensitized state of the channel and comprehensively describes the interactions of ligands of the orthosteric and benzodiazepines binding sites in agreement with the available experimental data. Furthermore, it is able to explain small differences regarding the binding of a variety of chemically divergent ligands. Finally, this new model may pave the way for the design of focused experimental studies that will allow a deeper description of the receptor.


Gamma-Aminobutyric
Acid type A receptors (GABA A Rs) are the main inhibitory neurotransmitter receptors in the mammalian central nervous system (Young & Chu, 1990). They are members of the Cys-Loop family of Pentameric Ligand Gated Ion Channels (PLGICs), along with the cation-selective, excitatory, nicotinic-acetylcholine receptors and serotonin receptors; and with the anion-selective, inhibitory, GABA C and Glycine receptors (Ortells & Lunt, 1995).
GABA A Rs are well characterized as pharmaceutical targets, and thus should be exhaustively studied: they have binding sites for a variety of ligands such as GABA, benzodiazepines, barbiturates, β-carbolines and neurosteroids among others (Werner Sieghart, 1995). In addition, they are involved in a myriad of neurological processes related not only to the regulation of inhibition in the CNS but also to the variability of GABAergic signals (Rudolph, Crestani, Hanns, & Rudolph, 2001;Vogt, 2015). Their correct functioning is extremely important for the health of humans; dysfunctional GABA A Rs have been related to anxiety, sleep disorders, epilepsy, alcohol dependence, among other affections (Collins et al., 2006;Crestani et al., 1999;Jones-Davis & Macdonald, 2003;Mukherjee, Das, Vaidyanathan, & Vasudevan, 2008;Nutt & Malizia, 2001;Yee et al., 2005) .
These receptors are generally divided into three distinct domains (Smith & Olsen, 1995): the extracellular, the transmembrane and the intracellular. The extracellular domain (ECD) is also known as the ligand binding domain because it hosts the orthosteric and benzodiazepines binding cavities. It is formed by 10 beta strands and two alpha helices. The transmembrane domain (TMD), which is formed by four alpha helices, controls the opening and closing of the channel pore through the movement of the inner (M2) helices. It also contains the binding cavities for other ligands, such as neurosteroids, ethanol and anaesthetics. The intracellular (IC) loop has not been completely solved experimentally, although it is known to be partially formed by an alpha helix. This domain helps modulate the flux of ions and the function of the channel by different mechanisms such as tyrosine phosphorylation (Kittler & Moss, 2003;Moss, Gorrie, Amato, & Smart, 1995) and the interaction with other proteins.
In this work, and due to their pharmaceutical relevance, we focused on the characterization of the orthosteric and high affinity benzodiazepines binding sites.
These cavities are homologous and they are located between different pairs of subunits. They are both formed by loops A, B and C from the principal subunits and loops D, E and F of the respective complementary subunits (Figure 1(a)). Different roles have been assigned to these loops in the binding of ligands and in allosteric interaction between sites, as well as in transmission of the signal in channel gating. Studies of the binding cavities of Cys-loop receptors and AChBPs with known structures (Calimet et al., 2013;Puthenkalam et al., 2016) have shown that they undergo conformational changes after the binding of the ligands. In particular, loop C has been shown to be very mobile, closing itself after agonist binding and adopting a more open conformation if the ligand is an antagonist (Puthenkalam et al., 2016). These cavities are allosterically and bi-directionally related: modulators that bind to the Benzodiazepines cavity modify the conformation of the orthosteric binding sites (Changeux & Edelstein, 1998) and it was also demonstrated that the binding of GABA and its activation of the receptor cause structural rearrangements in the benzodiazepines site (Teissére & Czajkowski, 2001). While the binding of agonists of the benzodiazepines binding site is positively coupled to the binding of GABA, imidazo-benzodiazepines are negatively coupled, i.e. these modulators stabilize different conformations of the receptor (Teissére & Czajkowski, 2001).
There are two orthosteric binding sites, which are located in the ECD between the principal face of a β-subunit and the complementary face of an α-subunit (β+/α-) (Figure 1(b)). While site G 1 is surrounded by a γand a β-subunit, site G 2 is flanked by an αand a γsubunit. Different compounds are known to bind to these cavities; those include the agonist GABA, muscimol, bicuculline and gabazine (SR-95531), (Figure 2). On the other hand, the allosteric high affinity benzodiazepines' binding pocket lies between the principal face of an αsubunit and the complementary face of the γ-subunit (α+/γ-) (Sigel & Buhr, 1997). The different compounds that bind to this site exhibit specific effects and binding affinities. Classic benzodiazepines (BZDs), such as diazepam, clonazepam, flurazepam and flunitrazepam, display a common 1,4-benzodiazepine nucleus with a 5phenyl substituent (Sternbach, 1979); imidazo-benzodiazepines (i-BZDs), such as flumazenil (Ro15-1788) and Ro15-4513, lack the 5-phenyl substituent but possess instead an imidazole ring substituted at positions 1 and 2 of the diazepine nucleus. Furthermore, there are other non-benzodiazepine ligands of this site, which include the imidazopyridine Zolpidem and the cyclopyrrolone Eszopiclone (S-zopiclone) ( Figure 3).
Acetylcholine-binding proteins and non-eukaryotic receptors such as ELIC (Erwinia ligand-gated ion channel), GLIC (Gloeobacter ligand-gated ion channel) and invertebrate Glutamate receptors have been extensively used in the past as templates to model mammalian Cys-Loop receptors (Berezhnoy, Gibbs, & Farb, 2009;Bergmann, Kongsbak, Sørensen, Sander, & Balle, 2013;Bertaccini, Yoluk, Lindahl, & Trudell, 2013;Carpenter, Lau, & Lightstone, 2012;Carpenter & Lightstone, 2016;Cromer, Morton, & Parker, 2002;Ernst, Brauchart, Boresch, & Sieghart, 2003;Hénin, Salari, Murlidaran, & Brannigan, 2014;Newell & Czajkowski, 2003;Puthenkalam et al., 2016;Sander et al., 2011;Thompson, Lester, & Lummis, 2010;Xie, Sha, Wang, & Cheng, 2013). However, in 2014 the first 3D structure of a GABA A receptor was published by Miller and Aricescu (Miller & Aricescu, 2015). The authors captured the spatial disposition of the atoms of a β 3 homopentamer, in a desensitized closed state, through X-ray diffraction with a 2.97Å resolution. The desensitized state is proposed to be of high affinity for ligands (Zhang, Xue, Liu, Yang, & Wang, 2013). This new structure is not only a priceless improvement for the modelling of these receptors, which in turn will provide new information to shed light on their structure and functionality, but also constitutes a prototype to analyse the binding of pharmacologically relevant compounds (Wang et al., 2017). Over the years, different techniques have been used to characterize the binding cavities, providing information about the location of the sites and the relevant residues involved directly in the interaction or affecting the efficacy of the different ligands. However, the lack of crystallographic structures of the α 1 β 2 γ 2 heteropentamer hampers the characterization of both the specific interactions in the cavity, and the mechanism underlying signal transmission and gating.
Here, using state of the art computational biology tools, we modelled the human α 1 β 2 γ 2 GABA A receptor. Subsequently, we characterized the interactions with the cognate ligands by performing docking experiments refined with molecular dynamics (MD) simulations. From the best of our knowledge, this is the first model able to widely explain experimental data on the aforementioned binding cavities. Indeed, it provides an extensive coverage of the relationship between experimental information and structural features of the orthosteric and benzodiazepines binding sites.

Constructing the model
The first step in the comparative modelling protocol consisted in retrieving the sequences of the α 1 , β 2 and γ 2 human subunits from Uniprot (Apweiler et al., 2004;Consortium, 2017): entries P14867, P47870 and P18507, respectively. For each subunit we searched for homologous sequences with SSearch (Pearson, 1991). The highest scores corresponded to two structures: the human GABA A Receptor β 3 homopentamer (PDB ID: 4COF) (Miller & Aricescu, 2015) and human Glycine Receptor α 3 homopentamer (PDB ID: 5CFB) (Huang, Chen, Michelsen, Schneider, & Shaffer, 2015). Although both of them were captured in the closed state, the former was co-crystallised with agonist benzamidine (desensitized state) while the latter was co-crystallised with an antagonist (closed-basal state). The differences between these structures involve the orientation of the M2 helices and the opening of loop C. The first receptor shows the highest percentage of identity with our sequences: approximately 40.8% with α 1 , 91.1% with β 2 and 43.2% with γ 2 , while the second shows 41.5% with α 1 , 45.4% with β 2 and 41.8% with γ 2 . We decided to carry out the modelling using the human β 3 GABA A Receptor as a template to prioritize a better quality of the sequence alignments. The percentages of identity are above the established threshold for confidence in the modelling (Table S3-SI) (Chothia & Lesk, 1986).
PROMALS (Pei & Grishin, 2007), HHPRED (Söding, Biegert, & Lupas, 2005) and Swiss Model (Arnold, Bordoli, Kopp, & Schwede, 2006;Biasini et al., 2014;Kiefer, Arnold, Künzli, Bordoli, & Schwede, 2009) were used to generate multiple alignments of the subunits with sequences belonging to other family members, and the crystallized homopentamer. The results from the three web-servers were compared to find differences only in loop F of α 1 subunit. Manual adjustments of two gaps in loop F of the α 1 subunits were applied in order to agree with experimental data (as from reference (Bergmann et al., 2013)).Truncation of intracellular loops was applied to the alignments following reference (Miller & Aricescu, 2015). The IC domain was not modelled due to lack of suitable templates. Indeed, this domain has not been established as involved in the binding of ligands in the ECD essentially due to the long distance between these domains. Naturally, it represents a limit in the implementation of the model to observe global conformational changes such as the opening and closing of the channel pore. 500 initial models were generated with Modeller 9.14 (Webb & Sali, 2002) using the 'Automodel class', which only includes spatial restraints obtained from the sequence alignment of the target and template, and required a refinement level 'refine.slow'. The best model was determined from the top 20 models, ranked according to the lowest value of the Modeller objective function and the 'Discrete Optimized Protein Energy' (DOPE) method score (Shen & Sali, 2006), as the one showing the largest percentage of residues in the most favoured region of the Ramachandran plot. ProSa z-score (Wiederstein & Sippl, 2007),Q-mean score (Benkert, Tosatto, & Schomburg, 2008) and PROCHECK (Laskowski, MacArthur, Moss, & Thornton, 1993) results were also taken into account in the selection process. These evaluation webservers have been extensively applied in several works concerning homology modelling (Carpenter et al., 2012;Roy & Mukherjee, 2017). Moreover, the top 20 models were carefully visually inspected to detect abnormalities. The chosen model was refined with Coot (Emsley, Lohkamp, Scott, & Cowtan, 2010) to optimize rotamers and side chain interactions, and PROPKA together with PDB2PQR (Dolinsky, Nielsen, McCammon, & Baker, 2004; webserver were employed to assign protonation states and optimize hydrogen bond networks.

Docking procedures
In order to find the binding modes of the selected ligands, we conducted a series of molecular docking simulations by two methods: blind and data-driven. We used these different techniques as a validation for the results when experimental information was not exhaustive enough.
The blind docking was performed with AutoDock Vina (ADV) (Trott & Olson, 2010), using a 25 Å x 23 Å x 20 Å grid centred at appropriate coordinates in each binding site. All the residues predicted by experimental data to be involved in the binding were within the gridbox. We prepared the ligands and the protein with Auto-DockTools (Morris et al., 2009) by adding missing hydrogens, combining the non-polar hydrogen atoms and computing Gasteiger charges. For each ligand a maximum of 20 binding modes were requested and the calculations were performed with the highest level of exhaustivity. The cut-off for the number of binding modes was a difference in the binding energy of 3 kcal/mol between the best and worst models. Both simulations with rigid and flexible side chains were performed, showing no considerable difference in the results.
On the other hand, we used HADDOCK 2.2, High Ambiguity Driven protein-protein Docking (Van Zundert et al., 2016), which employs experimental information about the interaction to establish the preferred binding modes. The experimental data, when available, was introduced as ambiguous interaction restraints (AIR constraints). The docking protocol implemented in the software is described elsewhere (Dominguez, Boelens, & Bonvin, 2003). The water refined structures were clustered using pairwise backbone Root Mean Square Deviation (RMSD) at the interface with a cut-off of 2Å, and analysed according to their average interaction energies (the sum of electrostatic, van der Waals and ambiguous restraints energies) and their average buried surface area. As stated in HADDOCK webserver, the parameterization of the ligands was done with PRODRG (Schüttelkopf & van Aalten, 2004).
We performed molecular docking simulations with four orthosteric ligands: two agonists (GABA and muscimol) and two antagonists (bicuculline and gabazine). The latter were studied although the receptor was modelled in an 'agonist-bound' state with the aim of assessing whether the model is yet adequate for the study of other classes of ligands. Following the same premise, we applied this technique to agonists, antagonists and inverse agonists of the benzodiazepines binding site; all of which are considerably larger than GABA or benzamidine.
The experimental information (Tables 4 and 5 -SI) used for the data-driven docking of each ligand was variable and strongly dependent on the availability. It should be stressed that, as implemented in HADDOCK, only a random 50% of these restraints were actually employed in the docking simulations. For residues of the orthosteric cavity we used the condition that they should interact with α 1 Arg67, β 2 Tyr97, β 2 Glu155, β 2 Ser156, β 2 Tyr157, β 2 Tyr205 and β 2 Arg207. However, no particular orientation was introduced as a restraint. As for the benzodiazepines binding cavity, the docking of Diazepam had very specific interactions incorporated as ambiguous restraints, specifically: the chlorine atom of Diazepam with α 1 His102 (Duncalfe, Carpenter, Smillie, Martin, & Dunn, 1996;Tan, Baur, Charon, Goeldner, & Sigel, 2009), and C3 of Diazepam with α 1 Thr207 and α 1 Ser206 of loop C (Tan et al., 2009). For the remaining ligands of the high affinity benzodiazepine's cavity only, the interacting residues were introduced as AIR constraints (Table 1).
The docking poses with the best scores according to both software programs were then analysed to compare the modelled interactions between ligand and receptor with those extracted from literature. The estimation of the energy of binding provided by the software programs was not contrasted with the experimental values due to a lack of direct correlation (refer to Table S8-SI for a discussion on this subject). Contacting atoms were defined as those separated by a distance shorter than the sum of their van der Waals radii plus 2.75 Å (approximate diameter of a water molecule). Afterwards, we identified through visual inspection those atoms that were engaged in biologically relevant interactions (i.e. through H-bonds, saltbridges, hydrophobic or cation-π interactions, among others). These interactions were tabulated and compared with those predicted by experimental studies, and from this comparison we computed two indices: Precision and Recall. These are statistical quantities that have been used in other works (Davis & Goadrich, 2006;Fierro, Suku, Alfonso-prieto, & Giorgetti, 2017;Goutte & Gaussier, 2005;Raghavan, Bollmann, & Jung, 1989;Saito & Rehmsmeier, 2015) to assess the quality of the binding modes resulting from docking experiments and to evaluate the performance of the techniques. They are, respectively, defined as: where TP 'True Positives' are residues predicted experimentally as relevant for the binding and found in our models to interact with the ligand; FN 'False Negatives' are the amino acids proposed by the literature to participate in the binding but are not interacting in our models; and FP 'False positives' are residues identified in our models as interacting with the ligands, although experimental data states otherwise, or are those amino acids within a contacting distance from the ligand but are not engaged in biological relevant interactions ( Figure S1-SI). The similarity of ligands might also shed light on their possible orientation inside the binding cavity, since it has been proposed that a relationship exists between structural similarity and a common binding mode (Boström, Hogner, & Schmitt, 2006). We employed OpenBa-belGUI, version 2.4.1 (Boyle et al., 2011;Morley, 2016) to calculate Tanimoto's coefficients to quantify the compounds' similarity ( Figure S4-SI).

Molecular dynamics simulations
Given the static nature of the docking protocols, the performance of molecular dynamics simulations would allow the ligands and receptor to adapt to the new conformation of the complex. It has been recently shown that in models of protein-ligand interactions of GPCRs, based on comparative modelling with low sequence identity, molecular dynamics simulations allowed for a better representation of the experimental binding configurations (Fierro et al., 2017;Gelis, Wolf, Hatt, Neuhaus, & Gerwert, 2012;Lai, Singer, & Crasto, 2005). We sequentially simulated the modelled receptor with different docked ligands: one molecule of GABA, muscimol, diazepam, clonazepam, flunitrazepam, flurazepam, Ro15-4513, flumazenil, eszopiclone and zolpidem.
The MD simulations were executed with GROMACS (Berendsen, van der Spoel, & van Drunen, 1995;Páll, Abraham, Kutzner, Hess, & Lindahl, 2015) using the SPC water model (Berendsen, Postma, van Gunsteren, & Hermans, 1981) and the GROMOS53A6 force field (Oostenbrink, Villa, Mark, & Van Gunsteren, 2004) for the protein. In all the simulations, the receptor, with the ligand attached, was embedded in a pre-equilibrated 1palmitoyl-2-oleoyl-phophatidylcholine (POPC) membrane and surrounded by a solution of water and CL−, NA+ ions in a concentration of 0.15 Mol. The Berger (Berger, Edholm, & Jähnig, 1997) parameters were used to characterize the lipids, and the parameters for the ligands were obtained from the Automated Topology Builder repository (Malde et al., 2011). The protocol consisted of an energy minimization of the system to a potential energy of 500 kJ/mol with Steepest Descent algorithm, secondly a thermalization in the NPT ensemble with the atom P8 of the phospholipids constrained to move on the plane of the membrane using the simulated annealing algorithm, then a protein/ligand-position restrained simulation to relax the interactions with the solvent, and finally the production runs. Van der Waals interactions were limited with a 1.2 nm cut-off and PME algorithm (Darden et al., 1993) was implemented for the electrostatic calculations also with 1.2 nm as cut-off. Verlet algorithm (Páll & Hess, 2013) was employed for neighbour search, and LINCS constraint algorithm (Hess, Bekker, Berendsen, & Fraaije, 1997) was used to restrain bonds' length. The temperature was kept at 310 K with Nose-Hoover thermostat (Hoover, 1985;Nosé, 1984) using a coupling constant of t T = 0.8 ps, and the pressure at 1 bar with Parrinello-Rahman barostat (Parrinello & Rahman, 1981) and using a coupling constant of t T = 10 ps and compressibility 4.5 10 5 bar. All the simulations consisted in production runs of 100 ns.

Model construction
The refined sequence alignment of each human subunit with β 3 -subunit, which was employed for the modelling process, is depicted in Figure 4. The extracellular and transmembrane domains consisted in 333 residues for α 1 and 332 for β 2 and γ 2.
The final model was chosen among the 500 initial configurations by considering structural parameters. The results of Q-mean programme (Benkert et al., 2008) showed that the model is slightly improved, compared to a previous model (Bergmann et al., 2013), and accept-
able relative to the crystallographic structure (Table S1-SI). In addition, PROCHECK scores were satisfactory (Table S2-SI) and the Ramachandran plot ( Figure 5) confirms the good quality of the backbone geometry, with 99.8% of the residues in the allowed regions.

Molecular docking
We performed blind and data-driven molecular docking for ligands of the orthosteric and high affinity benzodiazepines' binding cavity. The results were then statistically compared against experimental data available by calculating the recall and precision coefficients. The models that better agreed with experimental data were then funnelled through our analysis protocol.

Orthosteric binding cavity
We carried out docking simulations of four orthosteric ligands in both sites G 1 and G 2 (Figure 1(b) and S3). . Alignment of the human β 3 subunit (4COF) with human subunits α 1 , β 2 and γ 2 . The sequences numbering corresponds to the proteins in the mature form. The residues are coloured according to the percentage of identity with the consensus sequence. Dark blue represents an identity of 80% or more, blue corresponds to more than 60%, while light blue represents more than 40% identity. The loops that form the binding cavities are underlined in green, while the emblematic Cys-loop is underlined in yellow. The intracellular linker is underlined in light blue (ICL). The results presented here belong to site G 1 (except for the docking of bicuculline, which belongs to site G 2 ), since the observed binding modes were similar in both cavities. In general, all the studies showed very high precision and recall values (Table 2), indicating that there are few residues in the binding cavity that cannot be recognized by the docking protocols or are not available for interaction in our model. One exception is the binding of bicuculline, which shows a better recall value but lower precision. This might be caused by the size of this ligand, which occupies a large fraction of the binding cavity in disregard of the suggested interactions.
In the best docking mode obtained for GABA (Figure 6(a)), the ligand interacts through the amine-end with β 2 Glu155 and through the carboxylic group with α 1 Arg67, showing complete agreement with experimental information. Likewise, the result of the docking of muscimol ( Figure 6(b)) shows hydrogen bonds between the protonated amine with β 2 Glu155 and the backbone of β 2 Ser156 from loop B, and α 1 Arg67 with the ketone moiety. This mode agrees partially with previous models (Bergmann et al., 2013;Sander et al., 2011). We hypothesize that the difficulty in docking bicuculline might be due to the fact that this ligand is an antagonist, and therefore it is expected to bind to the cavity with loop C in a more open conformation, in contrast to our model, which represents an agonist-bound state of the receptor. In our best configuration (Figure 6(c)), α 1 Phe65 lies parallel to the benzene ring of bicuculline in a plane 3.6 Å below; they are possibly interacting via π-π stacking, in agreement with experimental data. Being smaller than bicuculline, gabazine is able to fit better inside the binding cavity (Figure 6(d)). Although the ligand interacts with α 1 Arg67 and β 2 Arg207 as expected, the experimentally suggested interaction with α 1 Asp63 is not present in our docking poses. Interestingly, bicuculline and gabazine displayed contacts with residues from the β-1 strand, namely α 1 Phe46, α 1 Val47, and α 1 Thr48, which, to our knowledge, have not been explored in experimental studies.
A complete list of the residues involved in interactions with the ligands of the orthosteric cavity in the best binding modes obtained through molecular docking is provided in Table 3.
The docking poses required further refinement since the model considered so far was only a static description of the receptor in a particular state. With this aim, we performed molecular dynamics simulations of the receptor in complex with GABA, as this is the most characterized agonist (Table 4 -SI), and also with muscimol. During the MD simulation, GABA relaxed in the binding cavity and formed (i) a cation-π interaction with an aromatic residue, i.e. β 2 Phe200, (ii) H-bonds with the tip of loop C, while maintaining the salt-bridges with β 2 Glu155 and α 1 Arg67 ( Figure S9-SI). The network of interactions with this compound was refined in such a way that precision remained in its maximum value, but the recall value improved (Table S6-SI). In addition, it can be appreciated that loop C opened slightly and loop F approached the principal subunit through residue α 1 Asp184, which came close to β 2 Arg207 (Figure S19-SI). In particular, residue α 1 Asp184 has drawn our attention, since it has been proposed as important in GABA binding, although, as it is located in loop F, it is far from the core of the binding site. However, this residue appeared in the docking results of Gabazine, the competitive antagonist, and it is observed that, while GABA remains bound, α 1 Asp184 moved towards the principal face of the cavity in molecular dynamics simulations.
β 2 Glu153 and β 2 Lys196 form a putative salt-bridge that regulates the movements of loop C during activation (Venkatachalan & Czajkowski, 2008). This interaction is supposed to be disrupted upon agonists binding and the consequent conformational change of loop C. In our initial model and throughout the MD simulation, although spatially close (~7 Å), these residues could not form a salt-bridge, which is consistent with the ligand-bound state of the template structure (Venkatachalan & Czajkowski, 2008). Moreover, in the cavity where GABA is bound, β 2 Glu153 is part of a stable interacting network which includes β 2 Glu155 and β 2 Arg207. On the other site, the MD simulation showed a change in the interactions, in which β 2 Glu155 now interacts with α 1 Arg120 and β 2 Arg207 ( Figures S20 and S21 -SI). This change in the interaction pattern might be related to the transmission of the agonist signalling.
The MD simulation of muscimol improved both scores to their maximum value (Table S6- β 2 Ser156, β 2 Tyr97, α 1 Thr130, β 2 Thr202 and α 1 Arg120 ( Figure S10-SI). The compound spent most of the simulation in a conformation similar to the starting pose.

Benzodiazepines binding cavity
Classical benzodiazepines share, among them, high similarity in their physico-chemical properties. These positive allosteric modulators are expected to bind near residues γ 2 Ala79 and α 1 His102, in particular the pharmacophore model derived from in vitro experiments (Clayton et al., 2007), suggested that the 7-substitutent of classical benzodiazepines should be located between residues Val203, Val212 and Tyr210 from loop C and His102 of loop A (Duncalfe et al., 1996) (Figure S5-SI). The best result for the docking of diazepam (Figure 7(a)) displays its chlorine atom pointing towards α 1 His102, α 1 Asn103 and γ 2 Asn60 and the C3 atom in the vicinity of the tip of loop C. In addition, the phenyl ring is surrounded by a hydrophobic cavity formed by α 1 Phe100, γ 2 Phe77, α 1 Tyr160, and α 1 Tyr210. As proposed in literature, the N-methyl substituent points outside the binding cavity (Tan et al., 2009) and the Note: Residues marked with an * are considered false positives. benzodiazepine ring lies above γ 2 Tyr58 in a parallel plane. Both precision and recall indices are high. Moreover, this configuration agrees well with those proposed by previous works (Berezhnoy et al., 2009;Middendorp et al., 2014) and with experimental information (Clayton et al., 2007;Table 5-SI). We compared this docking pose with those reported by Richter et al. (2012), finding our model resembles their second ranked mode, since it places the pending phenyl ring inside the lipophilic cavity in the back of the binding site, in contrast to their best binding pose which orients the ring outwards. The rest of the classic benzodiazepines also display high values of both precision and recall (Table 4). For these ligands, except Flunitrazepam (Figure 7(c)), the pending phenyl ring is surrounded by the characteristic hydrophobic cavity. Flunitrazepam's 7-substituent (a fluorine atom) interacts with α 1 His102 and γ 2 Tyr58, which cannot be completely fulfilled by flurazepam (Figure 7(d)), whose chlorine atom only interacts with γ 2 Tyr58 or by clonazepam (Figure 7(b)), which locates the nitro group above γ 2 Tyr58 and~7Å from α 1 His102. Remarkably, most of the interactions in the cavity are hydrophobic, which might account for the high values of the indices. Most of the residues involved in the binding of these ligands (Table 5) are the same as the ones suggested in the literature (Table S5-SI).
Bearing in mind the information available regarding the binding of diazepam, we performed two MD simulations, starting with the best configurations obtained from both docking methods. Throughout the simulations, Diazepam changed its orientation to adopt an almost identical configuration in the binding cavity for both MD simulations; the precision and recall scores exhibited small differences due to alternative disposition of the surrounding side chains. However, the final refined model showed improved indices (Tables S6 and S7-SI). On the other hand, along the MD simulations of clonazepam, this ligand sampled different conformations within the cavity so that the benzodiazepine ring lays parallel to the membrane plane, the nitro group points outside the cavity and the chloro-phenyl ring remains in the hydrophobic cavity. During the simulation with flunitrazepam the side-chains of surrounding residues relaxed in the cavity, while the ligand remained virtually motionless, improving slightly both evaluation parameters. The second best docking pose obtained for flurazepam (Figure 7(d)), although it shows lower precision and recall scores, resembles the pose adopted by bromoflurazepam while interacting with the prokaryotic receptor ELIC (Spurny et al., 2012) in one of its crystallographic structures (PDB ID: 4A98). Thus, we used this configuration to perform the MD simulations. Along the simulations, the ligand undergoes several changes in its orientation to finally converge to a conformation resembling the starting pose.
Regarding the imidazo-benzodiazepines (i-BZD), while they show high similarity between them, they are very different compared to classical benzodiazepines ( Figure S4-SI). Mutation and SCAM studies have shown that the 3' substituent of the i-BZD's imidazo-ring is located nearby γ 2 Ala79 and γ 2 Thr81 (Kucken, Teissére, Seffinga-Clark, Wagner, & Czajkowski, 2003), and the 7-substituent should be pointing towards α 1 His102, α 1 Val203, α 1 Val212 and α 1 Gly158 from loop B (Tan et al., 2007). α 1 Tyr210 is proposed to photoincorporate Ro15-4513, with the azide group pointing towards that region (Sawyer, Chiara, Olsen, & Cohen, 2002). The best docking configurations of Ro15-4513 (Figure 8.a) and flumazenil (Figure 8(b)) have their imidazo 3'-substituent pointing towards γ 2 Ala79 and the 7-substituent pointing towards α 1 His102. Additionally, α 1 Tyr210 interacts directly with Ro15-4513. Both models are in agreement with the requirements of the pharmarcophore model (Clayton et al., 2007). However, there are discrepancies regarding the orientation of the 4th and 5th positions of the benzodiazepine ring: these atoms point into the interior of the binding cavity in some of our models, but in one of the bests configurations obtained for Ro15-4513 the imidazole ring is orientated inwards, while these atoms point outwards. Previous models have oriented these atoms to the exterior of the cavity (Kucken et al., 2003) similar to the behaviour of the N-methyl and carbonyl moieties in Diazepam.
Since the binding modes for Ro15-4513 show identical scores (Table 4), we performed MD simulations for both systems. The model with the imidazole ring pointing inwards displayed better results: The simulation improved the precision and recall scores (Table 6), and the ligand adopted a stable configuration with loop C in an open conformation. Regarding the binding mode of flumazenil, MD simulations were not able to improve the docking precision and recall values. In our docking studies the best binding modes differed in the orientation of the imidazole ring. Considering that the MD simulations of Ro15-4513 gave better results for the system with the ligand oriented so that the imidazole ring points into the cavity, we believe that binding mode of flumazenil should also reproduce this feature.
The docking studies of zolpidem resulted in very similar conformations which displayed comparable recall and precision indexes ( Table 4). The observed interactions are mostly hydrophobic (Table 5). It has been already proposed that, since it does not share many polar interactions, the binding mode of zolpidem depends more on the shape of the cavity rather than on the interaction with particular residues (Hanson, Morlock, Satyshur, & Czajkowski, 2008). The best configuration (Figure 9(a)) can be related to previous computational studies performed with a different model (Vijayan, Bhattacharyya, & Ghoshal, 2012) and with the pharmacophore model based on experimental information (Clayton et al., 2007). However, in our model there is no direct interaction between the ligand and loop F as proposed in literature (Table 5 -SI). Table 5. List of residues that participate on the interaction with ligands of the high affinity benzodiazepines binding cavity of our modelled GABA A R.
Concerning the docking of eszopiclone (Figure 9.b), the best binding mode possesses high indices of precision and recall (Table 4) and agrees well with the pharmacophore model based on experimental data (Clayton et al., 2007). In addition, the orientation is similar to the one adopted by R-zopiclone bound to ELIC receptor (Spurny et al., 2012) (PDB ID: 4A97). Both ligands share common interactions with conserved aromatic residues (Table 5), even though the spatial disposition of the side-chains is not exactly pre-served. The chlorine atom is oriented towards the back of the cavity, between α 1 Tyr160 and α 1 Tyr210 and the methyl-piperazine group is located near α 1 His102 and α 1 Phe100, while the pyridine ring interacts with γ 2 Phe77 and α 1 Thr206. The pyrrolo-pyrazin is located under the tip of loop C and above γ 2 Tyr58. The free carboxyl interacts with γ 2 Tyr58 through and H-bond. A docking mode similar to that proposed by  was also found with ADV with the same binding score.   Both ligands remained bound during the MD simulations, with eszopiclone temporarily exploring other conformations along its simulation while zolpidem retained its orientation. Even though precision was increased for both ligands, recall was somewhat reduced (Table 6).
There are certain residues suggested in literature (Table 5-SI) as relevant in the binding of ligands but which are not observed in any of our models to have an active interaction with the compounds. For instance, residue α 1 Gly201 is located at the base of loop C; its mutation probably affects the binding of the ligands (zolpidem and eszopiclone) because longer chains might alter the movement capacity of this loop. This mutation might also modify the orientation of the neighbouring side-chains and directly affect the entrance of the ligand . Regarding the hypothetical interaction of flunitrazepam and zolpidem with γ 2 Met57, from our models we suggest that it may be very difficult for the side chain of this residue to interact directly with the ligands, because it is located in a beta strand with the side chain oriented outwards in the cavity. In addition, considering the following residue is γ 2 Tyr58, which has also been proved to directly interact with ligand, the disposition of the amino acids in the secondary structure seems correct. On the other hand, the experimentally observed reduction of binding affinities when α 1 Gly158 or α 1 Gly208 are mutated to Cys (Table 5-SI), can be explained by our model since a mutation in these positions would severely reduce the available volume of the cavity and alter the interactions between residues. The mutation α 1 G158C makes it possible for Ro15-4513 to interact with this residue through its azide moiety, as suggested experimentally (Tan et al., 2007).
Essentially, the refinement MD simulations performed on the receptor in complex with the allosteric modulators, improved both recall and precision indices (Tables S6 and S7-SI). All of the ligands remained bound to the cavity within the 100 ns of simulation and most of them did not change their orientation drastically (Figures S11-S18-SI). Flurazepam and diazepam displayed the most noticeable movements. The stability of the complexes and the ligands within the binding cavity was ascertained through measurements of the root-mean-square deviation (RMSD) (Figures S7 and S8-SI).
As a final remark, this work could provide a starting point for the design of several experimental studies. On the one hand, confirmation of the binding modes for the ligands is needed. On the other hand, we observed potential interactions between the ligands and unexplored regions of the receptor such as the amino acids from the β-1 strand in the orthosteric binding pocket, thus, we propose that mutational and substituted-cysteine accessibility method (SCAM) analysis would contribute greatly to understand the participation of these residues, which have been proposed to participate in the binding process in analogy with the benzodiazepines' binding site (Kucken et al., 2000). The approach of loop F of the orthosteric cavity towards the principal face through α 1 Asp184 in the MD simulations is a hint on the mechanisms involving this loop. Thus, we propose cross-linking and mutational studies should be performed on this residue in order to assess its role in the binding and activation of the channel, and the allosteric influence of benzodiazepines.

Conclusions
Due to their role in neurological health and vast pharmacology, the study of the structure of GABA A receptors, their function and interactions with other molecules is essential to understand the bases of many diseases and develop new treatments, as well as strategies for improving life quality.
It was not until 2014 that the first crystallographic structure of a GABA A receptor was published. For this reason, previous models were performed with templates sharing very low sequence identity belonging to other receptors of the Cys-loop family and AChBPs. In this work, we present a model of the α 1 β 2 γ 2 GABA A receptor in a closed-desensitized state based on the human β 3 homopentamer ( Figure S2-SI). The model here presented is able to reproduce active binding cavities, as they have been validated against all the existing experimental data.
Due to the lack of a heteromeric structure in the apostate or with bound ligands, state-of-the-art computational studies can be combined with experimental data to obtain information about the binding and conformational changes in the receptor that lead to its opening and closing. Docking results, although sometimes ambiguous, allowed us to propose binding modes for different ligands and helped in the formulation of hypotheses about the important residues. In this work, the docking procedures were able to capture the main features of the interaction of the receptor and its cognate ligands. Refinement using MD simulations allowed us to improve the results by providing a more extensive sampling of the conformational space of the ligand interacting with the receptor.
Regarding the binding properties of ligands of the orthosteric cavity, agonists GABA and muscimol could be docked in conformations that are similar for both compounds, in agreement with previous works (Bergmann et al., 2013;Boileau, Evers, Davis, & Czajkowski, 1999;Goldschen-Ohm, Wagner, & Jones, 2011;Holden & Czajkowski, 2002;Newell, McDevitt, & Czajkowski, 2004;Sander et al., 2011). Bicuculline, however, could not be located in a precise configuration in this model of the receptor, probably due to its size and the fact that it is an antagonist, and therefore it is expected to bind in the cavity with loop C in an open conformation. While gabazine, being smaller, fitted well into the cavity and was engaged in several of the suggested interactions (Boileau et al., 1999;Newell et al., 2004;Padgett, Hanek, Lester, Dougherty, & Lummis, 2007;Wagner & Czajkowski, 2001;Westh-Hansen et al., 1999) (Table 4 of the SI).
The current information about the high affinity benzodiazepines binding cavity allowed us to propose a binding mode for diazepam which is able to explain the available experimental data (Clayton et al., 2007;Middendorp et al., 2014;Tan et al., 2009;Table S5-SI). Other classical benzodiazepines, namely clonazepam, flunitrazepam and flurazepam, were docked in modes that are compatible with experimental data but which do not share exactly the same binding features as diazepam. The imidazo-benzodiazepines could also be docked based on experimental information (Kucken et al., 2003;Middendorp et al., 2014;Morlock & Czajkowski, 2011), showing that this site, in our model, is capable of binding not only positive modulators, but also antagonists and negative modulators. Non-benzodiazepine drugs zolpidem and eszopiclone exhibited binding modes that are also in agreement with experimental data . Remarkably, the docking of eszopiclone could be directly compared to the crystallographic structure of Rzopiclone bound to a prokariotic homologue receptor (ELIC) (Spurny et al., 2012). Both ligands show high similarity in the orientation within the cavity and make common key contacts with homologous hydrophobic residues of the binding sites.
In conclusion, we built a model of the α 1 β 2 γ 2 GABA A receptor, comprising both the extracellular and transmembrane domains, which was then docked with compounds of the orthosteric and benzodiazepines binding sites. All the described complexes are in agreement with the available experimental data. We, therefore, suggest that this model can be used in future studies related to the orthosteric and benzodiazepines binding sites. Additionally, we propose that mutational and SCAM analysis would contribute greatly to understand the participation of amino acids from the β-1 strand in the orthosteric cavity. Finally, we consider that this model might contribute to the future development of novel strategies for the design of specific orthosteric/allosteric ligands.