Exploring Volatile General Anesthetic Binding to a Closed Membrane-Bound Bacterial Voltage-Gated Sodium Channel via Computation

Despite the clinical ubiquity of anesthesia, the molecular basis of anesthetic action is poorly understood. Amongst the many molecular targets proposed to contribute to anesthetic effects, the voltage gated sodium channels (VGSCs) should also be considered relevant, as they have been shown to be sensitive to all general anesthetics tested thus far. However, binding sites for VGSCs have not been identified. Moreover, the mechanism of inhibition is still largely unknown. The recently reported atomic structures of several members of the bacterial VGSC family offer the opportunity to shed light on the mechanism of action of anesthetics on these important ion channels. To this end, we have performed a molecular dynamics “flooding” simulation on a membrane-bound structural model of the archetypal bacterial VGSC, NaChBac in a closed pore conformation. This computation allowed us to identify binding sites and access pathways for the commonly used volatile general anesthetic, isoflurane. Three sites have been characterized with binding affinities in a physiologically relevant range. Interestingly, one of the most favorable sites is in the pore of the channel, suggesting that the binding sites of local and general anesthetics may overlap. Surprisingly, even though the activation gate of the channel is closed, and therefore the pore and the aqueous compartment at the intracellular side are disconnected, we observe binding of isoflurane in the central cavity. Several sampled association and dissociation events in the central cavity provide consistent support to the hypothesis that the “fenestrations” present in the membrane-embedded region of the channel act as the long-hypothesized hydrophobic drug access pathway.


Introduction
Voltage gated sodium channels (VGSCs), which mediate the upstroke of the action potential in most excitable tissues, are key targets of anesthetics. The binding site and molecular mechanism of action for local anesthetics have been well characterized in the last few decades, while the role of VGSCs in general anesthetic action is less well understood and both mechanisms continue to be studied. Thus far, the VGSC binding sites for general anesthetics have not been identified. Identifying binding sites and access pathways for volatile general anesthetics is key to understanding their mechanism of action and to designing new drugs.
Sodium channels can be inhibited by a number of compounds, including toxins, quaternary ammonium compounds and local anesthetics [1]. When the pore is open, local anesthetics can enter from the intracellular side blocking conduction of ions. More recent work shows that general anesthetics also inhibit sodium channels, possibly through pore-blocking mechanisms [2][3][4][5][6][7]. Anesthetic block by both local and general anesthetics is preserved in the archetypal bacterial voltage-gated channel NaChBac [4,8].
Early observations [9] pointed to a local anesthetic binding site inside the channel, and indeed, in Nav 1.2, the binding site suggested by mutagenesis is in the central cavity [10]. The effect of local anesthetics on NaChBac is consistent with a cavity-binding site, although the precise binding pocket is probably not the same as in Nav 1.2 [8]. It is thus conceivable that local and inhaled general anesthetic sites in the channel cavity in NaChBac overlap.
The classical studies of charged local anesthetics and their analogs showed that blocking and unblocking seemed to require open channels. This led to the description of the ''hydrophilic pathway'' for drug access. However, additional experiments showed that hydrophobic local anesthetics could bind and unbind even when channels are closed [11][12][13]. This finding suggested an additional ''hydrophobic pathway'' that circumvents the closed gate. However, no structural correlate for this pathway has yet been conclusively identified.
Although large mammalian VGSCs have remained resistant to structural characterization, the discovery of the smaller bacterial VGSCs has provided a tool to characterize the structural features of these important channels [14,15]. In general, the major structural domains resemble those found in the crystal structures of voltage-gated K + channels ( Figure 1A). Transmembrane domains S1-S4 form the voltage-sensing domain (VSD), which is connected to the pore domain through the S4-S5 linker. The pore domain, formed by domains S5 and S6 includes the pore-loops (P-loops), which form the ion selectivity filter and may also constitute an inactivation gate [16,17]. On the intracellular side, four S6domains form the bundle crossing that constitutes the activation gate [18][19][20]. The S4-S5 linkers form a cuff around the bundle crossing, such that in the presence of a depolarizing stimulus, conformational changes in the VSD widen this cuff, allowing the bundle crossing to open to permit conduction.
A striking feature of the bacterial VGSC crystallized thus far is the presence of the so-called fenestrations [21][22][23]. These fenestrations provide a hydrophobic tunnel through the lipid-embedded portion of the channel to the central cavity, where the consensus local anesthetic site is located. In X-ray crystal structures, lipid tails occupy the fenestrations. These fenestrations may provide the hydrophobic drug access pathway that has been postulated for sodium channels since the 1970s [12,13,24].
The available atomic structure of a bacterial VGSC offers the first opportunity to address some of the most fundamental questions about the mechanism of volatile anesthetic action on VGSCs, namely what are the likely structural determinants of general anesthetic effects on VGSCs and could the fenestrations really provide access to the central cavity for small hydrophobic drugs? To investigate these questions, we present the results of a molecular dynamics (MD) simulation study of the bacterial VGSC NaChBac embedded in a lipid bilayer in presence of the inhaled general anesthetic isoflurane. Our unbiased ''flooding'' technique [25] is a way to flexibly dock isoflurane molecules to high affinity sites of the membrane-bound protein, as well as to thereby suggest energetically favorable drug binding pathways.

Identification of isoflurane binding sites
To identify putative isoflurane binding sites, we performed MD ''flooding'' simulations [25] on a structural model of NaChBac [26] inserted in a hydrated lipid bilayer (Figure 1, A and B, Movie S1). We considered a system comprising the membrane-bound channel along with a large number of drug molecules, initially located in the aqueous phase ( Figure 1B). As expected from the Meyer-Overton rule, we observed almost complete partitioning into the membrane within the first 100 ns of MD simulation ( Figure S1). In the subsequent MD trajectory, we observed binding to several hydrophobic pockets on the protein surface with at most a single isoflurane molecule in the aqueous phase. Importantly, the conformation of NaChBac was not significantly affected on the time-scale of our MD simulation by the binding of isoflurane molecules, as inferred from the RMSD of the channel along the MD trajectory ( Figure S2). The potentially large number of binding sites poses a challenge in identifying which of these may be pharmacologically relevant. To overcome this challenge we applied a two-fold strategy in which we first explored all the possible binding modes of isoflurane, and then used an established data-mining approach [27] to extract distinct configurations from our large data-set of molecular configurations. Specifically, we perform a cluster analysis on the positions of the centers of mass of the whole ensemble of isoflurane molecules and ranked the clustering solutions according to their local density. The rationale for this choice is to prioritize those regions that show continuous occupancy, which would suggest tighter binding of the drug molecule. This analysis revealed three major binding regions ( Figure 1C): 1) a region near the selectivity filter, termed the extracellular site, 2) a region near the S4-S5 linker, termed the linker site, and 3) a region within the cavity, termed the cavity site.
Due to the protein's four-fold symmetry, a functional channel has four equivalent copies of each of the three putative binding sites. To investigate whether the three putative sites are occupied symmetrically in the tetramer, and to determine which amino acids line each site, we analyzed the interactions between isoflurane and all the residues within the sites in each subunit. In particular, we monitored the contacts between any atom of the drug molecule and any atom of a given amino acid for each instantaneous configuration after isoflurane partitioning has reached equilibrium. Residues were classified as non-interacting, ''possibly'' interacting, and ''likely'' interacting according to the number of configurations in which the contact was detected. For the extracellular site, all amino acids in all four subunits are in contact with the drug. In the cavity site as well as through the fenestrations, most residues are possibly or likely interacting in all four subunits. Intriguingly, the linker site is asymmetrically occupied in two adjacent sites.

Characterization of isoflurane sites
We used the clustering data and proximity time analyses, in combination with analysis of hydrogen-bonding-like interactions and mobility of the drug molecule in the sites, to structurally characterize the sites and elucidate the key determinants of binding. Furthermore, to assess the pharmacological relevance of these sites, we estimated the binding free energy of isoflurane for each site using well-established FEP methods [28], allowing us to rank the relative free energy at each site. We also estimated binding affinities (Kd) and found that they were within a reasonable physiological range (see Methods) [29].
Extracellular site. The extracellular site sits at the intersubunit interface between the P-loops (

Author Summary
The molecular mechanisms mediating the pharmacologically induced state of general anesthesia are, in general, poorly understood. Modulation of voltage gated sodium channels is thought to play a major role in anesthesia, as several members of this class of channels show a significant response to general anesthetics. However, the detailed mechanism of inhibition or potentiation of these channels is completely unknown. Recently, the structures of several members of the bacterial family became available, thereby offering the opportunity to shed light on some of these issues. We have performed molecular dynamics simulations on one of these bacterial voltage gated sodium channels, NaChBac, to identify binding sites and access pathways for the volatile general anesthetic isoflurane. We found that isoflurane, at physiologically relevant concentrations, binds the channel at three distinct sites. One site is in the pore of the channel, suggesting that isoflurane may hinder the permeant sodium ions. Surprisingly, we found that this binding site is accessible to the drug even when the pore and the aqueous compartment at the intracellular side are disconnected. In our simulations, the ''fenestrations'' present in the membraneembedded region of the channel act as the longhypothesized hydrophobic drug access pathway.
duration of the simulation (Table S1). The estimated binding free energy (DG binding ) for this site is 24.260.8 kcal/mol. Consistent with the notion of a tight-binding interaction, rotational and translational diffusion in this site are restricted compared to aqueous solution (Table S2). Additionally, the corresponding diffusion constants are the smallest among all sites (Table S2).
Linker site. The linker site is located at the ''corner'' formed by the N-terminus of one linker and the C-terminus of the adjacent linker ( Figure 2, B and C). This site is the most dynamic and heterogeneous of the three sites, with variable occupancy and exchange of isoflurane molecules. Unlike the other sites, this site is asymmetrically occupied, with the a-d and c-d sites (Figure 2, legend) continuously occupied, while the other two equivalent sites are not occupied for any significant length of time over the course of the simulation (Table S1). Either one or two isoflurane molecules occupy the a-d site, while the c-d site is occupied by only one molecule. Also for this site, the majority of drug-protein interactions are non-polar (Leu133, Ala136 and Ile147) but a hydrogen-bonding-like interaction exists with Asn 225. The estimated DG binding for this site is 23.7620.4 kcal/mol for the first isoflurane molecule and 22.6620.9 for the second molecule, given the first one is already bound. In contrast to the extracellular binding site, only translational diffusion is significantly slower compared to aqueous solution, suggesting that isoflurane retains a significant degree of conformational flexibility inside the pocket.
Pore site. Though the activation-gate is closed, we find that the central cavity is occupied by up to five isoflurane molecules at once (Figure 2, D and E, Figure 3, A and B). All interactions in the pore are hydrophobic and isoflurane molecules associate with a number of residues. Importantly, we observed a stable interaction with residues known to play a role in local anesthetic action in mammalian sodium channels (F1764 and Y1771 in Nav1.2), which are partially conserved in NaChBac (T220 and F227) [8,10] (Table S1). Despite the fact that we observe exchange between the central cavity and the hydrophobic section of the lipid bilayer, the cavity is almost always occupied by at least one drug molecule (Figure 3 B). Therefore we expect the affinity for this site to be comparable to those for the other sites.
Given that isoflurane molecules in the central cavity are in direct contact with water, we were surprised that diffusion translational constant is comparable to those observed in the other two binding sites.
Fenestrations. Our homology model of NaChBac preserves the fenestrations found in the NavAb crystal structure [21,26], which provide a hydrophobic tunnel from the protein surface to the central cavity. We anticipated isoflurane occupation of hydrophobic pockets on the NaChBac protein surface, but were unsure whether isoflurane molecules would be able to fully traverse the fenestrations to displace bound waters and access the pore. In the flooding simulation, the S6 gate was closed, and there was no hydrophilic pathway for the drug molecules to the central cavity. Nevertheless, several isoflurane molecules were able to reach the pore through the fenestrations, displacing water molecules that left the central cavity through the selectivity filter. Interestingly, each isoflurane molecule is able to displace approximately six water molecules ( Figure 3A). Thus, in presence of a drug molecule, two solvation shells of water cannot surround a permeant sodium ion. An intriguing hypothesis raised by this observation is that inhibition may result from desolvation of the central cavity rather than steric hindrance, therefore a correlation between molecular volume of the drug and potency should be expected for this site. We found that isoflurane can both enter and exit the pore via fenestrations (Figure 4, Movie S2) and that all four fenestrations had both exchange and occupation by isoflurane (Table S1). This observation strongly suggests that fenestrations can act as a hydrophobic access pathway to the central cavity for small hydrophobic drugs such as isoflurane. However, it is still unclear if larger hydrophobic drugs, including other general anesthetics like sevoflurane and propofol, could traverse the same fenestration pathway. Crystal structures of NavBacs have so far shown that the fenestrations are occupied by the alkyl moiety of lipid molecules. Our simulations confirm that in the absence of isoflurane, each fenestration is occupied by the alkyl moiety of a lipid molecule; but in the presence of isoflurane, the lipid tails are displaced allowing isoflurane to enter the central cavity.

Discussion
Using an unbiased sampling of cavities within the bacterial channel NaChBac, we identified three putative binding sites for isoflurane. These binding sites are located at the extracellular side of the channel near the selectivity filter, at the linker between S4 and S5, and in the central cavity. Though the estimated binding free energies indicate high affinity, isoflurane molecules explore several orientations within each binding site. Due to the flexible nature of the Asn and Gln side-chains, isoflurane is able to form stable interactions with the protein at the extracellular and linker sites. Surprisingly, despite the fact that the pore is not accessible from the intracellular side, isoflurane molecules were able to reach the central cavity through the fenestrations. This observation suggests that the fenestrations could be the anesthetic hydrophobic pathways proposed more than three decades ago [12,13,24].

Possible mechanisms of anesthetic action
All three identified sites are in regions of the protein that are predicted by mutagenesis to be critical to gating and conduction, and we hypothesize that some or all of them will play a role in inhibition of NaChBac by isoflurane. Since the time-scales involved in the molecular events relevant for gating are beyond normal MD simulation time-scales, we cannot directly probe the effect of drug interactions on protein conformation. However, based on knowledge of the mechanisms of gating, conduction, and inhibition by local anesthetics, we postulate possible mechanisms of drug action involving the binding sites suggested by our simulations. Isoflurane binding to the extracellular site positioned in the P-loops could affect the conformation of the selectivity filter, leading to inactivation through filter collapse [16,17] and consequent reduction in peak current. Isoflurane binding to the linker site, which sits between the cuff formed by the S4-S5 linkers and the S6 bundle crossing, could affect the conformation or the coupling between the cuff and the bundle crossing to halt gating [30][31][32]. Finally, isoflurane in the pore site could simply occlude the pore through classical local anesthetic mechanisms [1], which has also been recently proposed as an inhibition mechanism of the GLIC channel by the same anesthetic [33]; this open channel block could be either purely steric or enhanced by the desolvation of the pore.
Though theories of anesthetic action focus primarily on anesthetic-protein interactions, the role of lipids cannot be discounted as partitioning into the bilayer changes lipid properties [34,35]. For the structurally unrelated channel gramicidin the effect of halothane on channel's function has been shown to be dependent on lipid composition [36]. Channel interactions with  lipids are important in regulating voltage sensing, gating and voltage-sensor-pore coupling [37][38][39] in voltage-gated ion channels. NaChBac's lipid sensitivity [40] suggests that appropriate lipid interactions regulate NaChBac function, and thus anesthetic partitioning could indirectly regulate the channel via disruption of lipid interactions. However, this simulation study is not designed to give realistic insights into lipid-anesthetic regulation for two reasons: first, because the lipid composition is not representative of a realistic bacterial or mammalian membrane and second, because the time-scale of the simulation is too short to observe lipid effects on channel structure and dynamics.

Testable hypotheses
This simulation of isoflurane binding sites and access pathways offers a number of experimentally testable hypotheses. While simulation shows that it is thermodynamically favorable for isoflurane to occupy these sites, it is possible that isoflurane occupation may have little or no physiological function. The functional relevance of each site can be ascertained by mutating each site individually to remove key determinants of isoflurane binding and evaluating whether NaChBac retains isoflurane sensitivity through electrophysiological assays [4]. The putative hydrogen bonding partners in the extracellular and linker sites (Gln 186, Asn 225 respectively) are particularly good candidates for mutagenesis, as point mutations changing the polar character of the residue are likely to cause a large change in binding affinity at these sites. Mutation of the cavity is likely to prove difficult, as isoflurane is highly mobile and interactions are nonspecific -the cavity site may be better probed through classical use-dependence and trapping assays [1]. Most intriguingly, the relevance of the fenestrations as a drug access pathway could also be probed through mutations that occlude the fenestration pore at the extracellular side.

Comparison to other anesthetic-channel simulations
Though no one has previously performed molecular dynamics simulation on voltage-gated sodium channels to identify anesthetic sites, numerous simulations have been done on ligand-gated ion channels (LGICs) due to the availability of the structures of bacterial homologues [41][42][43]. Simulation of anesthetic sites in LGICs, primarily based on the ELIC structure, have identified predominantly intersubunit sites, as well as a few intrasubunit sites [33,[44][45][46][47][48][49][50]. These sites are also the most common types identified experimentally [51,52], though pore sites have also been found [45,[52][53][54][55]. In fact, all isoflurane sites identified in NaChBac are intersubunit sites, with the pore site also being intersubunit. Importantly, the observation that binding occurs at the intersubunit regions in both voltage-and ligand-gated channels suggests that anesthetics may impair the cooperative conformational transitions that have been shown to be crucial for the function of ion channels.

Conclusions
Experimental verification of the aforementioned predictions will prove crucial to deepen our understanding of sodium channels modulation by anesthetics and will likely prompt further computational investigations. Indeed, despite providing a relatively accurate and detailed description of drug-channel binding events, our computational model is characterized by several limitations: (i) simulations were performed on a homology-based model of NaChBac; (ii) only one of the metastable structural conformations of the channel was probed for drug binding; (iii) lipid dependence of drug action was not addressed. Though the high degree of sequence identity with NavAb gives us confidence on the overall architecture of NaChBac, inaccuracies in the structure of the binding sites or major reorganizations of these pockets along the activation pathway can both potentially affect our results. Increased availability of experimental structures [22,56] and a detailed characterization of the conformational transitions entailed by the activation/deactivation cycle [57] will allow us to address these issues in future computational studies, as well as being able to begin to work with eukaryotic channels [58,59].

MD simulations
Simulations were initialized using the theoretical model of NaChBac in a closed conformation obtained previously [26]. This homology model was built on the basis of the X-ray crystal structure of NavAb (PDB ID: 3RVY), which was crystallized in a closed-pore conformation with all four voltage-sensing domains partially activated [57].
The model for the closed NaChBac channel, embedded in a fully hydrated lipid bilayer, was equilibrated by MD simulation. Specifically, the membrane is comprised of 1-palmytoyl-2-oleoylsn-glycero-3-phosphatidylcholine (POPC) lipid molecules. Isoflurane molecules were initially placed in the aqueous phase with random positions and orientation. The membrane protein complex contained a total of ,120,000 atoms, including NaChBac, 434 lipid molecules, 25310 water molecules, 236 ions in solution and 145 isoflurane molecules. The resulting initial aqueous concentration of isoflurane is 300 mM; the equilibrium aqueous concentration 0.9 mM. Two Na + ions were initially placed in the channel selectivity filter, in agreement with a previous computational study of NavAb showing double occupancy of the filter by Na + ions [60]. All charged amino acids were protonated using their respective pKas and the assumption that the solution was at pH 7. MD trajectories were collected for 0.5 ms.
MD simulation used the CHARMM22-CMAP force field with torsional cross-terms for the protein and CHARMM27 for the phospholipids [61,62]. A united-atom representation was adopted for the acyl chains of the POPC lipid molecules [63]. The water molecules were described using the TIP3P model [64,65]. Periodic boundary conditions were employed for all of the MD simulations and the electrostatic potential was evaluated using the particlemesh Ewald method [66]. The lengths of all bonds containing hydrogen were constrained with the SHAKE/RATTLE algorithm [67]. The system was maintained at a temperature of 300 K and pressure of 1 atm using the Langevin thermostat and barostat methods as implemented in the MD code NAMD2.8 [68] ( Figure  S1). The rRESPA multiple time step method was employed, with a high frequency timestep of 2.0 fs and a low frequency time step of 4.0 fs.
This computational setup has several limitations. First, the timescales achievable by our MD simulation are too short to observe channel gating. Second, the small size of the lipid bilayer in the simulation box results in an excessive drug concentration in the bilayer. The final limitation stems from the use of a oversimplified POPC lipid bilayer. Bacterial sodium channel function is strongly lipid dependent [40]. However, NaChBac is functional in liposomes containing POPC and in mammalian cells, making POPC a reasonable choice for NaChBac simulations.

Cluster analysis
To identify binding sites, we analyzed equilibrated configurations of the system and seeking regions characterized by high density of isoflurane. We first computed the center of mass (COM) of each isoflurane molecule in the MD trajectory frame. We then performed a cluster analysis on the resulting set of COM positions using the geometric distance between each pair of positions to build a proximity matrix. Partitioning of the set (clustering) was obtained using the Jarvis-Patrick algorithm [27] with a nearestneighbor list of 8 and shared-neighbor threshold of 3. We then ranked the clusters according to their average density. Here, we treated each COM position as an isotropic Gaussian density of width 1 Å , and summed over all the COM positions belonging to a given cluster. Integration of the density for each cluster was performed on a 3-dimensional grid with bin dimensions of 0.560.560.5 Å 3 . The top three ranked clusters, comprising 15% of the total number of configurations, are discussed as putative binding regions.

Free energy calculations
Free energy calculations were performed using the free energy perturbation (FEP) method. The binding free energies are calculated using the following scheme: DG bind = DG gas-prot 2 (DG solv +DG rstr ). Here, DG bind is the free energy of binding isoflurane to NaChBac, DG vac-.prot is the free energy of transferring an isoflurane from the gas phase to the binding site, DG solv is the isoflurane solvation free energy, and DG rstr is a measure of the entropy cost associated with the reduction in volume from a 1 M solution (V 1M ) to the volume available at the binding site (V rstr ), i.e., DG rstr = RT ln(V rstr /V 1M ). For all reported binding energies, V 1M is given by the volume associated with the flat-bottom spherical restraint applied to keep the isoflurane in the binding site during the interaction decoupling. Calculations were performed in NAMD 2.8 by varying the coupling parameter in steps of 0.025 at the ends and 0.05 in the middle.
This approach has been successfully applied to binding of anesthetics to proteins including the binding of R-and Sisoflurane enantiomers to apoferritin [69] as well as the binding of isoflurane and propofol to the GLIC bacterial ion channel [33]. The estimated binding free energy for isoflurane was assuming that the only two relevant thermodynamic states are proteinbound isoflurane aqueous isoflurane. Therefore the water-lipid partitioning is not considered. However, since the affinity of isoflurane for a hydrophobic phase is significantly lower than our estimated affinities (DG of transfer from water to dodecane is 3.0 kcal/mol), we expect the lipid partitioning to have a marginal effect on the apparent K d . Movie S1 Isoflurane flooding. The NaChBac pore domain (blue ribbons) sits in the bilayer (headgroups shown as grey spheres) while isoflurane (red/green molecules) begins in the aqueous compartment above and below the bilayer and partitions first into the bilayer and then into the protein structure. Voltage sensing domains, water molecules and lipid alkyl tails are not shown in this movie but are present in the simulation ( Figure 1A).

(MPG)
Movie S2 Representative isoflurane trajectory through fenestration to cavity site. Bottom view of NaChBac structure (grey ribbons) surrounded by lipid molecules (stick representations).
Isoflurane (blue, space-filling representation) enters from the lipid phase into the cavity through a fenestration. (MPG) Table S1 Proximity of isoflurane to residues within each site. The table reports the relative frequencies (as a percentage of frames) of contacts between isoflurane and the residues lining each binding site (a contact is detected if the distance between any atom of the isoflurane molecule and any atom of the given residue is smaller than 5 Å ). Residues are classified as non-interacting (red), possibly interacting (orange), and likely interacting (green) depending on the probability of being engaged in a contact. (DOCX)