A Review of the Updated Pharmacophore for the Alpha 5 GABA(A) Benzodiazepine Receptor Model

An updated model of the GABA(A) benzodiazepine receptor pharmacophore of the α5-BzR/GABA(A) subtype has been constructed prompted by the synthesis of subtype selective ligands in light of the recent developments in both ligand synthesis, behavioral studies, and molecular modeling studies of the binding site itself. A number of BzR/GABA(A) α5 subtype selective compounds were synthesized, notably α5-subtype selective inverse agonist PWZ-029 (1) which is active in enhancing cognition in both rodents and primates. In addition, a chiral positive allosteric modulator (PAM), SH-053-2′F-R-CH3 (2), has been shown to reverse the deleterious effects in the MAM-model of schizophrenia as well as alleviate constriction in airway smooth muscle. Presented here is an updated model of the pharmacophore for α5β2γ2 Bz/GABA(A) receptors, including a rendering of PWZ-029 docked within the α5-binding pocket showing specific interactions of the molecule with the receptor. Differences in the included volume as compared to α1β2γ2, α2β2γ2, and α3β2γ2 will be illustrated for clarity. These new models enhance the ability to understand structural characteristics of ligands which act as agonists, antagonists, or inverse agonists at the Bz BS of GABA(A) receptors.

. Absolute subunit arrangement of the  GABA(A) receptor when viewed from the synaptic cleft. 25 The GABA(A) binding sites are located at the     subunit interfaces and the modulatory benzodiazepine binding site is located at the    subunit interface.
The 2 subdomain is important for interacting the benzodiazepines in addition to many other molecules that act via the Bz BS. 11,12 Sequence variations that differ between  and γ subdomains are responsible for the subtype selectivity and efficacy of Bz BS ligands. [11][12][13][14] Three segments have been proposed to generate the Bz BS: 1) The   side, 2) the loops A, B and C, and 3) three segments of the   called loops D, E and F. 11, 12 X-ray crystallography and EMstructures AChBP 1 and the nAChR 15 confirmed these segments, which create a groove-like pocket between subdomains and that are conserved among the entire superfamily. Trudell and Hagen proposed a near planar cleft BzR binding sites in 1980. 16,17 Figure 2. Alignment and homology model depiction of the so called "loops A-F" and flanking regions of the human sequences of the Bz recognition site in different subunits. 1,[11][12][13][14][15]25 View of the homology model ( Figure 2) with the channel in between units. The membrane can be pictured parallel to the top edge of the structure. Stick representation is used for key residues: Loop A His102, loop B Tyr 160 and loop D Phe77.
After the nAChR structure was available, newer reports of liganded AChBP crystal structures were published. 18 Importantly, these structures showed that the local conformation of the BS is influenced by the ligand. This is particularly the case for the loop C. Thus, this highly mobile structure element responds towards ligand binding and induces smaller changes along its boundary. 18 This phenomena has been observed for other receptors as well and it has been hypothesized that low energy conformers can be separately stabilized by different ligands. Unfortunately, a priori prediction of the appropriate conformer/ligand pair is not conclusive at this time. The models can predict changes in binding sites as large as 40% and distances of key residues can vary several Ångstroms.
Small changes of the protein conformation can severely change the efficacy of Bz BS ligands. 19 In addition, GABA(A) receptor ligands have different efficacies towards particular GABA(A) receptor subtypes. 20,21 The stabilization of the active state can be as small as 1 kcal/mol underlining the fact that proteins are inherently dynamic and can adopt many conformations. 22 Therefore, the prediction of the absolute assignments of specific side chains to specific descriptors for any particular conformation is not possible at present. However, a "suggested" orientation of the pharmacophore in the receptor can be predicted implying conformational flexibility of residues that satisfy pharmacophoric descriptors. The combination of a pharmacophore model and a homology results in a model that assigns large lipophilic areas to specific regions and enables flexible assignment of for instance H-bond bridges or π-π stacking. Water molecules can significantly complicate prediction of binding. Therefore, our group favors the pharmacophore/receptor models for drug design based on the rigid planar ligands (K i = 5nM-20nM) that have introduced by Trudell et al. 17

Relative Orientation of the Pharmacophore within the Comparative Model
The Milwaukee group published a review prior to structure determination of the AChBP that included results of site-directed mutagenesis and gave valuable insights in respect to important side chains that formed the Bz BS and represent pharmacophoric descriptors. 23 The combination of that model and the recent homology models to lymnea AChBP, 2 aplysia AChBP and the nAChR 24 resulted in a novel model that is discussed here. Experimental evidence induces some degree of conformational flexibility of the model that hint towards variable assignments for H-bridge interactions. In addition, specific areas of lipophilic interactions create a binding site geometry flanked by particular amino acid residues ( Figure 3).

Combining Homology Model and Experimental Evidence
Docking of diazepam and flumazenil was performed to test our computational ligand docking to homology models as described by the two assignments above. Different models based on different templates were used 1,15,18   All poses were found in the database. Representatively, diazepam is shown in Figure 5.  25 In Figure 5, diazepam was shown as space filling structure (turquoise: carbon, white: hydrogen, green: halogen, red: oxygen, and blue: nitrogen). The receptor was presented at ribbon with key amino acids in stick representation. The insert figure represents the empty pocket in the "upright" position. This orientation is different from the main figure, where the structure has been turned and tilted to bring diazepam into the orientation of the unified pharmacophore model. In this orientation, His102 and Thr207 satisfy H 2 and H 1 , respectively, and His 102, Val 212 and Tyr 210 would represent L 2 a hydrophobic pocket. This orientation of the 7-substituted reactive compounds such as isothiocyanate (see Figure 4, third structure) is in position to react with His102, Tyr 210, Val 212, and Val 203. This orientation was developed in collaboration with Professor Werner Seighart. 25 Figure 6. High resolution image of diazepam docking pose in binding pocket. 25 The docking pose in Figure 6 and similar predictions satisfy the two definitions for L 2 and A  that are discussed above and the following descriptors: H 1 is aligned with the side chains near the top of loop C, the L DI region is close the subdomain interaction (the 1 loop B (Y160) and the  region are forming the sheet involving M130, T142 and F77); A 2 is represented by Hbond interaction near 1Y160, A79, M130 and T142. If these proximity relations that are found in the docking poses are compared with the unified pharmacophore model we can predict the following orientation the residues.  25 The Figures 7 and 8 show the orientation of Bz BS residues relative to the pharmacophore model in two views that correspond to a 90° rotation representing the evaluation of experimental data and the comparative model of the  GABA(A). The ligand is diazepam. It is possible that more than one side chain could satisfy the same descriptor because inverse agonists stabilize protein conformations that vary from the conformations preferred by agonists, In addition, there is still some degree of variation in structural details of the poses that represent this orientation. The factors of these structure variations: template and alignment choices, construction of the model and its refinement, and docking algorithms.
With the alignment of XLi-093 and docked Ro15-1788, 25 it is becomes clear that these compounds are able to orientate towards the extracellular domain instead of solvent accessible space through the subunit boundary. Orientations of the pharmacophore model inside the structural model of the GABA(A) receptor have been proposed in the past. [28][29][30] Current refinements of the predicted protein-ligand complexes with the support of the unified pharmacophore model is a promising approach to create 3D structures that are more accurate and applicable for structure-guided drug design. New approaches that combine protein modeling and pharmacophore elucidation are currently investigated.

Construction of the Unified Pharmacophore/Receptor Model
The 150+ BZ ligands 31, 32 that belong to fifteen different scaffolds were used to create a unified pharmacophore/receptor model. The relative affinities, efficacies and functional effects from ligands of the same structural class at the diazepam-sensitive and diazepam-insensitive benzodiazepine receptor binding sites applied. In addition, we took the approximate locations of descriptors (hydrogen bond donor sites, hydrogen bond acceptor sites, lipophilic regions, and regions of steric repulsion) that were based primarily on in vitro binding affinities in account. The compounds from different scaffolds were superimposed in order to satisfy the same descriptors, creating the unified pharmacophore model this included the rigid planar templates (Ki = 5nM-20nM) of Trudell. 17 The model has two hydrogen bond donating groups (H 1 and H 2 ), one hydrogen bond accepting group (A 2 ) and one lipophilic site (L 1 ). There are also lipophilic regions of interaction (L 2 , L 3 , L 4 , and L Di ) and regions of negative steric repulsion (S 1 , S 2 and S 3 ). Important for positive allosteric modulation are occupation of L 2 and/or L 3 , and L 4 as well as interactions at H 1 , H 2, and L 1 . Inverse agonists only require interactions with the H 1 , L 1, and A 2 descriptors for good activity in vivo. [33][34][35][36][37] The L 3 region is lipophilic of nature differentiating diazepam sensitive (DS) and the diazepam insensitive (DI) coordination. Figure 8 depicts the relative locations of the different descriptors and regions and classes of ligands for GABA(A) R are shown in Table 1.   . Alignments of several Bz BS ligands within the pharmacophore model. 25 A unified pharmacophore model incorporating many substance classes that act at the DS and DI benzodiazepine binding sites of GABA(A) receptors has been updated to include new substance classes. Compound development guided by this pharmacophore model has led to new agents with interesting pharmacological profiles, particularly enhanced preference for  or  containing GABA(A) receptor subtypes. Based on the evaluation of experimental data and comparative models of the  GABA(A) receptor, the location of several residues relative to the descriptors of the pharmacophore/receptor model has been proposed. Although no absolute assignments were made regarding which amino acids satisfy the pharmacophoric descriptors, experimental data strongly indicated definite trends with regard to how ligands of varying pharmacological activity are oriented within the receptor. Because the unified pharmacophore/receptor model accounts for the binding and activity profiles at the six GABA(A) receptor subtypes containing any one of the different alpha subunits, the proposed orientation should also be similar within the different models 1 of the various receptor subtypes. Information to be immediately gained from these proposed orientations can have far reaching benefits, not only for the rational design of selective ligands and the interpretation of ligand docking results, but also for the identification and evaluation of possible roles certain residues may have within the pocket. As structure determination of the GABA(A) receptor is eagerly awaited, it is hoped that these proposed orientations may be used by others to gain additional insight into the potential mechanisms underlying binding and modulation at the Bz site, all of which will lead to a better understanding of the structure and function of GABA(A) receptors. In the next section the methods used to build the protein model will be explained.

Homology Models of the Benzodiazepine Receptor
The -aminobutyric-acid (GABA(A) and GABA C ) pentameric ligand gated ion channels are members of a superfamily of allosteric transmembrane proteins which includes the nicotinic acetylcholine (nAChR), serotonin 5-HT3, and glycine receptors. Electrophysiological data on GABA(A) is available but attempts at atomic resolution to acquire structural data have so far been unsuccessful. However, in 2001 Smit 38 discovered a homologue of the nACh receptor ligand-binding domain from the snail Lymnaea stagnalis. Acetylcholine-binding protein (AChBP) is produced in the glial cells of Lymnaea stagnalis. In the synaptic cleft, AChBP modulates synaptic transmission. Subunits of the snail glial cells form a stable homopentamer with conserved N terminal domains. Known agonists and antagonists of nAChRs bind to the homologue. For this reason, it is a valuable template for modeling the N-terminus domains of pentameric ligand gated ion channels.
Brejc et al. successfully solved the crystal structure of AChBP by X-ray crystallography using weak Pb multiple wavelength anomalous diffraction (MAD) data in two crystal forms. 1 The crystals were grown at room temperature using the hanging-drop vapour diffusion technique. The crystal structure revealed a radially symmetric homopentamer with extracellular dimensions of the nicotinic acetylcholine receptor (nAChR) which correlated with measurements taken using electron microscopy. 1 Figure 10. Orthogonal views of the homopentameric acetycholine binding protein crystallized from Lymnaea stagnalis. The subunits have been assigned different colors for clarity. 2 Subunits are homologous and have been colored by chain to distinguish them. In space filling models, HEPES (N-2-hydroxyethylpeperazine-N'-2-ethanesulphonic acid (present during the crystallization process) can be seen bound in the ligand binding site ( Figure 10) between each subunit interface. 2 The N terminus is located at the "mouth" of the ion channel in the synaptic cleft. A molecule of HEPES was present in the acetylcholine binding pocket. This was verified as the residues implicated with agonist binding in nAChR and AChBP were conserved in the receptor. The residues lie in four regions identified as loops A, B, C, and D. Loops A, B, and C are located on what is referred to as the 'plus' face while loops D, E, and F lie on the 'minus' face (see Figures 11 and 12). In the heteromeric nAChR, agonist binding takes place between the 'plus' face of a subunit and the 'minus' face of an  or γ subunit. Site directed mutagenesis and radioligand binding assays have established that residues of the GABA(A) receptor involved in ligand binding align with residues in loops A, B, and C of the  subunit and loops D and E of the subunit. 2, 23, 25  Brookhaven National Laboratories, this is the home of the Protein Data Bank. The PDB archive is the repository of information about the 3D structures of large biological molecules such as proteins and nucleic acids. With the protein sequences in hand, a critical first step in building a homology model is to align the raw protein sequence with the protein sequence of a target protein or template structure. Alignments are prepared using multiple sequence alignment software for DNA and proteins. 39 After calculating the best match for the selected sequences, they are lined up so identities, similarities, and differences can be seen. Programs which perform this include MAFFT, Muscle, Multalin, ClustalW, and BlastAlign (http://pbil.univ-lyon1.fr/alignment.html). Templates can also be found by performing a query of the raw protein sequence against the Protein databank. This is performed using the protein database search program, Gapped Blast. 39, 40 The report will contain a list of sequences producing significant alignments, PDB code, protein description, a normalized alignment score and an Expect (E) value. In general, a lower E score describes the background noise that exists between sequences. A lower E value correlates to a more significant score. Once a proper alignment of protein sequences was performed, the raw protein sequence was threaded onto the AChBP using Deep View. The crystal structure of the AChBP was downloaded from the Swiss Protein Databank Viewer. Deep View is a Swiss Protein Databank Viewer (http://ca.expasy.org/spdbv). Threading is performed after the alignment is complete. Amino acids which have been identified as equivalent between the proteins are "fit". This is similar, but more precise than the "3 corresponding atoms" technique. Deep View includes a tool, Magic Fit, 39 which performs the threading procedure for all residues across a raw sequence and a reference protein. Once the preliminary model has been built, a manual inspection can be performed. Considerations during manual inspection include, identifying areas of low homology, and manually aligning gap regions with loops in the reference protein. 39,40 Loop regions which correspond to gaps were in the alignment can be modeled by fitting structures from a loop database. In the current model, loop regions were mapped to a loop database by Cromer et al. 2,15 Upon completion of the preliminary model, the following properties can be studied using Deep View and Sybyl 25 : 1. Surface hydrophobicity 2. Solvent accessibility of N-glycosylation sites after pentamer assembly 3. Acidic,basic, and non polar maps 4. Sybyl X minimizations 5. Steric clashes 6. RMSD of alpha carbons 7. Location of disulfide bonds In order to improve the prediction of areas of low homology PHD was utilized by Cromer et al. 2,39,40 PHD or more commonly referred to as PredictProtein is a sequence analysis tool and the prediction of protein structure and function. Protein sequences or alignments are submitted to Columbia University; PredictProtein returns multiple sequence alignments, PROSITE sequence motifs, low-complexity regions (SEG), nuclear localisation signals, regions lacking regular structure (NORS) and predictions of secondary structure, solvent accessibility, globular regions, transmembrane helices, coiled-coil regions, structural switch regions, disulfide-bonds, sub-cellular localization, and functional annotations. Upon request the Rost group will perform fold recognition by prediction-based threading, CHOP domain assignments, predictions of transmembrane strands and inter-residue contacts are also available. PredictProtein is run by Burkhard Rost at Columbia University in New York (www.rostlab.org). 40 As an alternative, Swiss Model can be used to model subunit loop regions by searching structures from loop databases.
Once subunits were completed, a pentamer was constructed taking into consideration residues determined experimentally to be present in the allosteric and GABA binding site. On the  subunit, this includes Tyr160, Tyr210, Phe77, and Ser205. 23,25 After manually assembling the subunits and comparing an overlay with AChBP, the structure was energy minimized. 25,39 With protein receptor models for and completed a truly unified pharmacophore model was assembled for each receptor subtype. Using the Tripos Biopolymer software, the pharmacophore was manually aligned into each protein. With mutagenesis data on residue interaction on key ligands, the pharmacophore was carefully adjusted so as to recreate the protein-ligand interactions discussed previously. The protein models below represent the cumulative results of over 20 years of research. The proteins are oriented with the  on the left and the  subunit on the right with loop C reaching over the middle to . Diazepam has been docked in each receptor and the pharmacophore is visible.

Protein Models of the  and Subtypes
Presented here are the complete unified protein-pharmacophore models. The Milwaukeebased pharmacophore has been inserted into each of the homology models for the receptor subtypes of the benzodiazepine receptor. The models show the similarities we expected. The two oxygens of serine 206 represent the hydrogen bond acceptor A 2 . The hydrogen bond donor, H 1 , is due to the hydroxyl group oxygens of Threonine 231 and Tyrosine 234. Threonine 193 is responsible for the H 2 hydrogen bond donor as its hydroxyl group coordinates with lone pairs of ligands. The imidazole ring of histidine 102 coordinates with phenyl rings of ligands through  stacking interactions.  From this unified model which was constructed on binding data from hundreds of compounds, included volume analysis, site directed mutagenesis, manual as well as Sybyl Flexidock and AutoDock algorithms, one can now see how any compound will dock in the benzodiazepine binding site of GABA(A) . Presented below is a compound for further research on drug abuse, WYS8.

Figure 28. Flexidock fit of WYS8
Refinement of WYS8 in the unified pharmacophore receptor was completed using FlexiDock. The protein with the pharmacophore set inside the benzodiazepine receptor is displayed with a ligand docked to the pharmacophore. From the compute function, Flexidock is executed. The pocket needs to be defined for Sybyl if it has not already. Pick atoms which are located near the binding site. All atoms within 4 Angstroms will be automatically selected by Sybyl to assist. For Flexidock to work with rigid proteins the following steps must be taken: 1. Water must be removed 2. Select atoms around the binding pocket 3. Add hydrogens 4. Pick 4Å radius to display around binding pocket 5. Have Sybyl add charges if necessary 6. Set rotatable bonds From these models, the individual pharmacophores, and the compound database the medicinal chemist has a powerful set of tools to determine future ligands to synthesize in order to deliver subtype selective drugs with reduced side effects. For example, knowing how to dock ligands into the pharmacophore, one can begin taking measurements of ligands and residues and analyzing interactions. In the Bz receptor His 102 and the phenyl ring of diazepam have a  stacking distance of 4.1 Angstroms. Next we will discuss how one can properly dock ligands into the binding site.

Docking Studies Using AutoDock 4.2
Docking of ligands into the new binding sites of the homology models was also executed using the suite of Autodock Tools developed at Scripps Institute. 41, 42 The general procedure entails preparing pdb file formats of the receptor (protein) as well as the ligand. Once these docking grids and docking parameters are set, the docking run will produce multiple potential docking reports. Prior knowledge from SAR studies and mutagenesis studies is then used to select those for a parallel run (clusters) and a final docking prediction is developed.

Preparing the Ligand for Autodock
The ligand is loaded into Autodock as a pdb file. Once in Autodock, Gasteiger partial charges are determined along with aromatic carbons, rotatable bonds and torsional degrees of freedom. A new "pdgqt" file is created which has this information attached to it. Preparation of the receptor begins with loading it too as a pdb file into Autodock. Water molecules must be removed from the protein if present for Autodock processing. If any residues are flexible, these need to be identified in the receptor. Once identified, torsion angles can be assigned. Bonds in side chains can be manually or automatically selected as rotatable, unrotatable, or non-rotatable. Autodock needs two files when docking ligands in proteins with flexible residues, a flexible PDBQT file and a Rigid PDBQT file. Before docking can begin the receptor is checked for Gasteiger charges. A map is created which identifies possible hydrogen bonding interactions.
Next, a grid box ( Figure 29) is constructed around the general vicinity of the receptor if this is known.

Figure 29. Setting the autogrid in Autodock for the Bz BS protein
With all the appropriate information stored in a parameter file, the docking program will know which map file to use, details about the ligand rules, the potential presence of flexible residues and which docking algorithm to use. Four different docking algorithms are currently available in AutoDock. Monte Carlo simulated annealing (SA) was used. Briefly, the ligand makes random moves in the receptor and the energy of the new position is compared to the old one. Moves in which the energy is lowered are accepted. WYS8 was drawn in ChemDraw Ultra 8.0 ( Figure 30). Lone pairs were added and the molecule was minimized in ChemDraw 3D using MM2 to a minimum RMS gradient of 0.100. This was converted to a PDB file to be loaded into AutoDock 4. Upon loading into Autodock, Gasteiger charges were added, 14 non-polar hydrogens were merged, 10 aromatic carbons were identified, 3 rotatable bonds were detected, and TORSDOF was set to 3.

Figure 31. Rotatable bonds are identified along with root atom in ligand WYS8
The root atom with the smallest largest sub-tree (in this case an acetylenic hydrogen) is marked with a green sphere and the root portion (red) is identified. The rotatable bonds are colored green. The receptor was rigid for this analysis. Using Tripos Benchware 3D Explorer 2.6, one of the results was modeled ( Figure 32). Key residues were identified and the protein surface five angstroms from the ligand was mapped showing its lipophilic potential. The results using AutoDock and a rigid protein will be much improved when flexible residues are allowed in the binding pocket along with rotatable bonds. The docking in Figure 28 is preferred when considering the orientation in the Milwaukee-based pharmacophore.
Homology models which have the pharmacophore aligned in the binding pocket have been presented for the and  GABA(A) receptors. The binding site was constructed in light of years of work including structure-activity-relationships (SAR), site directed mutagenesis, and docking studies of monomer and dimers alike. The and receptors have not been made as these "diazepam-insensitive" binding sites are not a focus of our research to date. An updated included volume analysis has been presented using ligands from 15 different structural families. The new included volume analyses and the unified pharmacophore protein models will continue to drive the design of novel and selective benzodiazepine receptor ligands. From these new ligands, scientists will increase their understanding of the physiological responses of the specific receptors, the design of pharmaceuticals with reduced side effects, and treatments for schizophrenia, drug addiction and Alzheimer's disease. The receptors are of particular interest because of their concentration in the hippocampus. The discovery of a larger L 3 pocket in the  receptor versus the other subtypes will provide medicinal chemists with a structural advantage in designing more subtype selective ligands. From the unified pharmacophore protein models presented here, the design of new compounds for the GABA(A) will be more productive than ever.

Future Work
Screening and modeling of more subtype selective ligands will further our current understanding of the benzodiazepine receptor and its function. Continued work on the molecular spreadsheet which will include all the compounds designed in the Milwaukee laboratory as well as other organizations will create a powerful QSAR algorithm for answering the age old question, "What compound should we make next?" This study will need to include all binding and non-binding compounds in order to create a predictive program with a high degree of accuracy (r 2 , q 2 ) in predicting activity.
Site directed mutagenesis needs to be continued in which receptors of all subtypes are tested to determine how a residues substitution affects each of the subtypes. Such experimentation will shed further light on how ligands are positioned in the binding pocket and what residues they are interacting with. Adding to the new pharmacophore protein alignments which can now be executed, docking runs using Autodock in which the residues in the receptor are set to flexible may give further information as to docking preferences of ligands helping to better understand key interactions in each of the different subtypes.
All of these activities have allowed us to develop the current ligands as well as propose next generation ligands. The new pharmacophores will aid in furthering our understanding of GABA and in designing better drugs for the CNS until a complete GABA(A) protein receptor ion channel is captured in high resolution.   1 1 1 0 1 1 1 1 0 1 0 0 1