Assessment of Homology Templates and an Anesthetic Binding Site within the γ-Aminobutyric Acid Receptor

FOR more than 160 yr, anesthetics have been used safely and effectively to minimize the otherwise deleterious side effects of major invasive procedures. Despite such success, their exact mechanism remains elusive. This is confounded by the fact that our understanding of the conscious state which general anesthesia alters is grossly inadequate. However, over the past years, great insights have been gained into the molecular underpinnings of anesthetic mechanisms, just as they have with comparable factors within the central nervous system. Anesthetics are thought to mediate a significant portion of their activity via binding and modulating the transmembrane ligand-gated ion channels (LGICs). In particular, the γ-aminobutyric acid receptor type A (GABAaR) and the glycine α1 receptor (GlyRa1) are the ion channels whose inhibitory currents are potentiated by the presence of general anesthetics. It is the GABAaR whose perturbation by barbiturates and benzodiazepines produces a state strongly similar to general anesthesia. Likewise, the intravenous general anesthetics, propofol and etomidate, are thought to modulate consciousness through specific sites in these channels.¹  These ion channels are composed of five subunits, often heteropentameric and variable in stoichiometry, arranged around a central ion-conducting pore. The extracellular domain of a representative subunit is characterized by a large component of β-sheet secondary structure and contains the binding site for the native ligand germane to the channel in question (i.e., γ-aminobutyric acid for the GABAaR). The transmembrane domain is composed largely of four-helix bundles. The latter secondary structure is a significant feature predicted by our group and later validated by experimentalists.^2,3  It is within this transmembrane region that anesthetics and alcohols are thought to bind and convey the majority of their channel modulation.⁴ 

However, because the LGICs are transmembrane proteins, their isolation and purification have proven difficult. To date, high-resolution crystal structures of the GABAaR do not exist. However, significant progress has been made toward understanding their structure and interactions with anesthetics via the methods of molecular modeling. The mainstay of such calculations lies in the techniques of homology modeling. Homology modeling is the method by which the amino acid sequence of a protein of unknown structure is aligned and threaded over that of a closely related amino acid sequence with known three-dimensional structure, such that the coordinates of the known protein can be transferred to those of the unknown. While such homology modeling involves a great deal of computational theory, it is also very dependent on experimentally described coordinates of proteins to act as templates with high sequence homology to the desired protein.^4,5  Over the last several years, several templates with great homology to the LGICs have been determined via cryo-electron microscopy, x-ray crystallography, and nuclear magnetic resonance that have made model construction more robust. In this article, we described how our new models based on such a template can account for much of the currently available experimental data concerning these channels, allowing correlation of ligand-binding measures with experimental potencies. These models have been used to illustrate the mechanism of channel gating^6,7  and can now possibly be used as the bases for future in silico high-throughput screening for new anesthetic discovery.

The cryo-electron microscopy–derived structure of the α1 subunit from the torpedo nicotinic acetylcholine receptor (nAChR, PDB ID 2BG9),⁸  the nuclear magnetic resonance–derived transmembrane domain structures of the α4 and β2 subunits of the nAChR (PDB ID 2LLY, 2LM2),⁹  and the crystallographically derived structures of the eukaryotic glutamate-gated chloride channel (GluCl, PDB ID 3RIF, 3RHW),¹⁰  as well as the prokaryotic pH-sensing channels, Gloeobacter violaceus ion channel (GLIC, PDB ID 3EAM),^11,12  and Erwinia chrysanthemi ion channel (PDB ID 2VL0),¹³  were downloaded from the Research Collaboratory for Structural Biology protein Data Bank. Most of the calculations below were performed with algorithms present within the Discovery Studio 3.1 software suite (Accelrys, San Diego, CA) except where noted. Multiple structure alignments were separately performed via the SAlign¹⁴  and 3DMA algorithms. A multiple sequence alignment profile was then derived from these structural alignments. The sequences of the individual subunits for the mouse GABAaR were downloaded from the National Center for Biotechnology Institute. They were concatenated and arranged in the clockwise order (as viewed from the extracellular side of the channel) γ 2, α1, β3, α1, β3 in accordance with our previous work¹⁵  and the experimental work of others.^16–19  These works demonstrate the appropriate stoichiometry and arrangement of subunit subtypes so as to satisfy a variety of mutational results required for specific ligand binding in the extracellular domain of the GABAaR.^20,21  The sequences of the mouse GABAaR α1, β3, and γ2 subunits were then aligned with the previously derived sequence of the known structural template (GluCl) using the Align123 algorithm, a version of the ClustalW routine²²  (see Supplemental Digital Content 1, https://links.lww.com/ALN/A955, which is a figure containing the color-coded multiple sequence alignment showing the individual subunits from the templates, target sequences, and several other homologous sequences within the LGIC with annotation for homologous amino acid positions relevant to anesthetic modulation).

The GluCl structural coordinates were used as the template for homology modeling for the reasons reviewed in the Discussion section. For the GluCl template, the native glutamate ligand, which was located in the intersubunit cleft within the extracellular domain, and the ivermectin located in the intersubunit cleft in the transmembrane domain were removed. The Modeller²³  algorithm was used for the assignment of coordinates for aligned amino acids. Modeller implements comparative protein structure modeling by satisfaction of spatial restraints that are obtained from the set of structural templates, thereby generating a probability density function governing the way in which heavy atom coordinates are assigned. Routines implemented within the Modeller framework also allow the construction and optimization of loops of amino acids, for which no template coordinates can be assigned, and the initial refinement of amino acid sidechains. Final rotamer search optimization was performed using the ChiRotor routine within Discovery Studio.

A second GABAaR model was independently constructed using comparative modeling in the Rosetta software suite.²⁴  The new eukaryotic GluCl structure is likely the one structure closest to GABAaR, but this channel has an ivermectin molecule bound between subunits. To work around this, we used the GluCl structure for modeling of each subunit, and the prokaryotic GLIC channel for the relative placement of the subunits. Each GABAaR subunit sequence (α, β, γ) was aligned separately to GluCl, with the long M3-M4 loop in the GABAaR subunits replaced with six glycine residues. The GluCl structure (PDB ID 3RHW) was used as a template to create 10,000 models with fragment-based modeling, and then clustering was used to select a model for each subunit. These subunits were superimposed on the GLIC structure in clockwise order γ, α, β, α, β (as viewed from the extracellular side), after which the sidechains where repacked using the SCWRL4 algorithm.²⁵ 

Due to the great similarity between the two independently built GABAaR models as signified by the root-mean-square deviation analyses noted in the Results section, a single model (from the Modeller routine) was used for all subsequent docking computations. As in previous work,⁴  the binding sites were readily identified by the convergence of residues relevant to anesthetic modulation on a common transmembrane intersubunit location in three-dimensional space. The volume of the initial binding site was very small due to the aforementioned energy optimizations and the rotation of amino acid sidechains into the binding site. Therefore, a single propofol molecule was manually inserted into one of the two putative binding sites between α and β subunits within the model built using Modeller. This propofol was optimized into position using molecular mechanics with the CHARMm force field, followed by molecular dynamics,²⁶  while allowing for sidechain flexibility with only backbone atoms fixed. Although the work by Roberts et al. has demonstrated that occupancy of only one binding site is sufficient to enhance the homologous glycine receptor’s function,²⁷  the work by Forman et al.^28–31  shows that propofol’s enhancing action in GABAaR requires occupancy of perhaps three propofol sites. For the sake of computational simplicity and because both binding sites are identical in amino acid contribution and location between α and β subunits, only one transmembrane binding site was chosen for docking examinations. The CDocker algorithm³²  was then used to randomly dock a congeneric series of flexible propofol derivatives, as well as propofol itself, into a rigid binding pocket and score their best relative binding affinities for correlation with known GABAaR potentiation EC₅₀s. CDocker is an implementation of a CHARMm molecular mechanics-based docking tool using a rigid receptor. This method begins with generating a set of 10 ligand conformations using 1,000 steps of high-temperature (1,000 K) molecular dynamics with different random seeds and includes electrostatic interactions. Random-pose orientations of the conformations were produced by translating the center of the ligand to a specified location within a predefined sphere incorporating the receptor site and performing a series of random rotations. The sphere in this case was chosen as a 7-Å radius from the center of mass of the initial manually docked propofol pose. This was done in an effort to both center the region of docking exploration within the area of relevant residues and to have a large enough radius to allow adequate random placement of these residues. A softened (temporarily reduced van der Waals radii) energy was calculated and the orientation kept if the energy was less than a specified threshold (300 kcal/mol). This process continued until either the desired number of low-energy orientations was found (set at 10) or the maximum number of bad orientations were tried (set at 800). Each orientation was subjected to simulated annealing molecular dynamics. The temperature was heated to a high enough level (700 K over 2,000 steps) to overcome local energy barriers between adjacent minima and then cooled to the target temperature (300 K over 5,000 steps). A final minimization of the ligand in the rigid receptor using a nonsoftened full potential (including full van der Waals radii and electrostatics) was performed. For each final pose, the CHARMm energy (interaction energy plus ligand strain) and the interaction energy alone were calculated. The poses were sorted by CHARMm energy and the top scoring (most negative, thus favorable to binding) poses were retained. For performance, many of these steps used a nonbonded energy grid, rather than the full potential energy terms usually used by CHARMm. This choice can provide a significant time saving during large database screening, but was reverted to its fully explicit CHARMm form during final optimizations to avoid the loss of some initial accuracy associated with the initial grid-based search.

The compounds with agonist activity that were docked included propofol (2,6-diisopropylphenol); 2,6-dimethylphenol; phenol; 2,6-diethylphenol; 2-isopropylphenol; and 2,6-disecbutylphenol. This is the same series of propofol analogs studied by Eckenhoff et al. for their correlation of binding free energy to horse spleen apoferritin derived from isothermal titration calorimetry with experimentally derived EC₅₀ for GABAaR potentiation.^33,34  These compounds were specifically chosen for our docking studies because they are the only propofol derivatives which also have experimentally derived specific binding constants to horse spleen apoferritin. It is the latter which were initially and reliably reproduced with the CDocker methodologies as a means of method validation specific to this set of compounds (unpublished data). Finally, as inactive controls, both the cis and trans forms of the nonimmobilizer 1,2-dichlorohexafluorocyclobutane (F6) were docked to the same binding site. This compound is without known potency for response to a painful stimulus (hence the name, nonimmobilizer) and is experimentally devoid of activity at the GABAaR.³⁵  As additional controls, we also docked the following propofol analogs noted for their inactivity at the GABAaR in the study by Krasowski et al.³⁴  such as dicyclohexylphenol, ditert-2,6-butylphenol, cyclopentylphenol, and dicyclopentylphenol. Linear regression analysis correlating the best negative CDocker score of each compound with the negative logarithm of its experimentally derived GABAaR EC₅₀ for potentiation was performed within Microsoft Excel 2007 (Microsoft Inc., Redmond, WA).

The sets of structural alignments from both the SAlign and 3DMA algorithms (see Supplemental Digital Content 1, https://links.lww.com/ALN/A955, and Supplemental Digital Content 2, https://links.lww.com/ALN/A956, which are figures containing the color-coded multiple sequence alignments) revealed that there was clear consensus as to the spatial location of residues relevant to anesthetic effects based on all templates, except for the intermediate resolution structure of nAChR (2BG9). That is, when multiple sequence alignments are performed between all of the structures noted previously along with the GABAaR and the GlyRa1, the multiple structure alignment agrees with the positions of homologous residues except for the structure of 2BG9. These issues with the 2BG9 structure have been noted by others as well.^10,36,37  For the anesthetic effect in the human GlyRa1, these positions are those that are homologous to residues ILE 229, SER 267, and ALA 288. The 2BG9 coordinates had gross discrepancies of amino acid positions of homologous residues known to alter anesthetic modulation in several transmembrane domains. The specific misaligned residues were dependent on the structural alignment algorithm, but in neither case could the location of all relevant residues be completely satisfied by the 2BG9 structure. In particular, using the SAlign algorithm, in the α1 subunit of the nAChR, it was LEU 257 on TM2, and LEU 279 on TM3 that did not align with their homologous positions on the remaining templates, despite multiple sequence alignment to the contrary (figs. 1 and 2). Optimizing to a slightly different structural alignment with the 2BG9 template using the 3DMA algorithm, it was ILE 220 on TM1 and LEU 257 on TM2 that could not align with their homologous positions on the remaining templates. This result was the justification for not using 2BG9 as a template for subsequent modeling, in addition to the further rationale for using GluCl noted in the Discussion section.

Once the GABAaR models were constructed, the preponderance of secondary structure in the extracellular native ligand-binding domain remained in the β-sheet conformation, whereas that of the transmembrane domain remained as tetrameric α helical bundles. This was noted as a consensus in the models derived using both the Modeler and the Rosetta methods. In fact, the Modeller routines were run in the facilities in Palo Alto, whereas the Rosetta routines were run in Stockholm, Sweden. Despite their independent constructions by completely different methods and in different locales, the models proved extraordinarily similar across the entire pentameric assembly as well as in the important putative transmembrane anesthetic binding clefts. In particular, the root-mean-square deviation of the α and β dimer, between which anesthetic binding is thought to occur, for the two models is on the order of 1.56 Å based on main chain atomic alignment using the SAlign routine. Of particular importance, the root-mean-square deviation values for the three critical residues involved in a transmembrane anesthetic binding site (β N265, β M286, α L231) are very comparable as follows: 0.597 Å for α carbon atoms, 0.599 Å for all main chain atoms, 1.426 Å for all atoms. As noted in our previous work for the GlyRa1, the majority of residues in homologous proteins that were experimentally labeled as facing hydrophilic regions also lined the water accessible portions of our model structures, whereas those which were experimentally labeled by lipophilic reagents were often accessible to the lipid membrane portion of our model.⁴  Likewise, experimental data on cysteine crosslinking of adjacent residues may be partially accounted for with these models as they are consistent with the ability to form disulfide bonds between specific intersubunit residue positions on the αTM1 subunit and the adjacent βTM3 subunit.^38,39 

The three residues in the mouse GABAaR, which were at homologous positions and notable for modulating anesthetic action within the transmembrane domain (β3 N265, MET286, and αL231), continue to line the intersubunit interface between α and β subunits after energy optimization (fig. 3), as noted in previous models.⁴  Again, this similarity was noted as a consensus in the models that were independently derived using both the Modeller algorithm and the Rosetta methods. Docking of the propofol derivative series to the region bound by these residues showed strong linear correlation with experimentally derived GABAaR potentiation EC₅₀ (fig. 4) as suggested by an R² of 0.92.

CDocker scoring was also validated by the results of the control compounds with known inactivity at the GABAaR. The CDocker scores for dicyclohexylphenol, ditert-2,6-butylphenol, and cyclopentylphenol were −14.3, −13.2, and −6.0 kcal/mol, respectively. These latter three are all significantly worse than the lowest potency agonist studied (phenol with a CDocker score of −16.5 kcal/mol and an EC₅₀ of 920 μM), suggesting the latter three predicted potencies are well into the millimolar range and beyond the ranges experimentally tested by Kraswoski et al.⁴⁰  The Cdocker-binding score was very unfavorable for dicyclopentylphenol (6.5 kcal/mol). The so-called “F6” nonimmobilizer from the Eger group was docked in both its cis and trans conformations. The CDocker-binding energies were grossly unfavorable for both versions of this compound (cis = 29.1 kcal/mol and trans = 28.4 kcal/mol) to the putative anesthetic binding site. This is in stark contrast to the CDocker energy scores for all of the propofol derivatives studied here that are agonists and which are clearly in the energetically favorable negative range.

For many years, homology modeling of the GABAaR has been hampered by a limited amount of highly homologous structural templates. In 1985, the first low–resolution cryo-electron micrograph of the homologous nAChR was released.⁴¹  Over the next 20 yr, this model was refined to the intermediate resolution of 4 Å and published in the Research Collaboratory for Structural Biology Protein Data Bank as entry 2BG9.⁸  This was a tremendous accomplishment at the time. However, although the 4-Å cryo-electron micrograph provided great insight into the overall secondary, tertiary, and quaternary structure of the class of LGICs, the exact atomic delineation of individual amino acid assignments in three-dimensional space had been called into question. After all, this cryo-electron micrograph was the consensus of micrographs over the four different subunit subtypes within the torpedo nAChR (two α, one β, one γ, and one δ). At 4-Å resolution, the globular appearance of sidechain identities is insufficient to make precise identifications of individual amino acids. Although this was the best technology could offer at this time, it is still amazing how much correct information is, in fact, present in this structure.

Given the recently released nuclear magnetic resonance structures of the transmembrane domains of the α4 and β2 nAChR subunits and their high homology to the torpedo nAChR, it is fairly clear where the discrepancies exist in the amino acid assignments within the subunits of torpedo nAChR (2BG9), making this a much less desirable modeling template from which to build a GABAaR structure.⁴² 

The GluCl and GLIC templates were initially considered for the homology modeling based on their homology to the GABAaR sequence, the availability of a whole pentameric ion channel construct, and their low root-mean-square deviation from the nuclear magnetic resonance structures of the α and β subunits of the nAChR. The use of templates only with all five subunits (i.e., GLIC and GluCl) arranged in a pentameric symmetry is vital for the analysis of anesthetic interactions because the location of this critical anesthetic binding site within the GABAaR is actually between specific α and β subunits within the transmembrane domain. Furthermore, the GluCl structure is thought to be in the open state of the ion channel pore, a feature important for anesthetic binding, because the anesthetic potentiation of the GABAaR is believed to come from the binding to and stabilization of the open channel state.²⁹  It is also the GluCl structure noted in the 3RIF and 3RHW templates that actually shows a ligand, ivermectin, bound in the same transmembrane location as an anesthetic binding site in the homologous GABAaR. Ivermectin is known to produce a similar potentiation for the GluCl as anesthetics do for the GABAaR. This GluCl structure also has the native bound ligand, glutamate, which further tends the construct toward an open structure. In addition, this structure would be similar to the condition of anesthetic potentiation of the GABAaR in which its native ligand, γ-aminobutyric acid, must be present. Furthermore, GluCl is the only template derived from the eukaryote, Caenorhabditis elegans, as opposed to GLIC, which is derived from the prokaryote, G. violaceus. Because eukaryotes have cholesterol in their membrane structure and prokaryotes do not, GluCl may be more representative of the proper intersubunit helical packing found in higher eukaryotic organisms. For these reasons in the end, the GluCl template formed the basis for modeling in this article.

When residues that are known to modulate anesthetic activity from homologous positions in different LGICs are mapped to the sequences of known template structures, consensus structural alignment based on the remaining five templates excluding 2BG9 revealed the consistent formation of an intersubunit anesthetic binding cavity within the transmembrane domain of these LGICs (fig. 3). This consensus lends greater confidence to the overall model building process and to further pursuit of validation analyses. Such analyses have included the ability to simulate the large-scale motion of the channel that was performed in our previous efforts in which we described an “iris-like” motion of channel gating.^6,43,44  Further validation also came in the form of illustrating channel responses to ligand binding via long-term molecular dynamics simulations. This has been accomplished with molecular dynamics simulations in our previous models of the GlyRa1 in the presence and absence of ethanol for 1 μs periods.⁵  However, actual channel gating occurs on timescales that are orders of magnitude larger and will require computational resources not currently available. Additional validation is now presented in this article via the correlation of relative binding energy scores with experimentally derived measures of GABAaR potentiation.

In the past, there has been modest controversy regarding the location of anesthetic binding sites within the LGICs. In this article, our goal was to examine one of those sites. There is evidence that anesthetics bind to sites within the transmembrane subunits of homologous proteins, specifically in GLIC and the E. chrysanthemi ion channel. In former models of various LGICs based on older templates, certain constructs were presented which postulated such an intrasubunit-binding site for anesthetics within the glycine receptor as well.⁴⁵  However, newer more homologous templates allow the construction of more robust models in which there is a convergence of residues relevant to anesthetic modulation on a common intersubunit binding site. As suggested in our recently published article,⁴⁶  a key site for anesthetic and alcohol action may be between the transmembrane subunits. An even more interesting occurrence has been shown in the crystal structure of E. chrysanthemi ion channel with the anesthetic bromoform in which both potentiation and inhibition have been noted by binding of the anesthetic to an intersubunit transmembrane cleft⁴⁷  in addition to effects within the pore and extracellular domains. Furthermore, the extracellular domain within this class of ion channels involves its own set of inter- and intrasubunit clefts. Pan et al. have demonstrated that the intravenous general anesthetic, ketamine, binds at an intersubunit site and causes inhibition of the channel in the extracellular or native ligand-binding domain of the protein. However, our intention in this work is to merely examine one of these putative binding sites which has been extensively studied in our previous modeling works and the experimental work of other groups.^4,46,48,49 

Once the binding sites can be identified as such, further validation and use of models can also come from predictions of anesthetic affinity and possible activity in previously known compounds, which have heretofore not been known to have anesthetic potency. This goal would be accomplished through the application of a wide variety of molecular docking and scoring schemes which allow individual compounds to be docked to a designated anesthetic binding site within the protein in question, as has been carried out by Liu et al.⁵⁰  Any such endeavor would be based on the ability to at least predict the orders of potency within a series of similar anesthetic ligands from a qualitative standpoint and, when feasible, predict a quantitative measure as well. The energetic parameter which should ultimately correlate with such potency is the free energy of ligand binding. However, total free energy is dependent on both enthalpy and entropy. Accounting for the binding entropy changes can be very difficult using computational methods. Ligands which vary markedly in size and shape could have very different entropic contributions to binding that are poorly accounted for by anything other than very complex and computationally expensive free energy perturbation calculations. Such calculations are impractical for high-throughput screens. However, within a reasonably similar congeneric series of ligands, as used here, such an entropic contribution can be assumed to be similar amongst members, with the enthalpic contribution of binding dominating the relative interactions. It is this enthalpic contribution that is most closely related to the docking scores presented here. With this, our model of an anesthetic binding pocket has characteristics that allow reasonable correlation of ligand enthalpic docking scores with actual experimentally derived measures of GABAaR potentiation for the set of propofol derivatives studied (fig. 4). Although our models will certainly undergo further refinements as new data become available, the current construct is among the first models to be validated by such a correlation. Although molecular mechanics may predict binding reasonably well, it remains to be seen whether all binding results in correlated efficacy/agonism. This process may allow the identification of possible new anesthetic ligands via methods of in silico high-throughput docking screens and perhaps be of benefit to future anesthetic drug design.