Molecular Simulations of Hydrophobic Gating of Pentameric Ligand Gated Ion Channels: Insights into Water and Ions.

Ion channels are proteins which form gated nanopores in biological membranes. Many channels exhibit hydrophobic gating, whereby functional closure of a pore occurs by local dewetting. The pentameric ligand gated ion channels (pLGICs) provide a biologically important example of hydrophobic gating. Molecular simulation studies comparing additive vs polarizable models indicate predictions of hydrophobic gating are robust to the model employed. However, polarizable models suggest favorable interactions of hydrophobic pore-lining regions with chloride ions, of relevance to both synthetic carriers and channel proteins. Electrowetting of a closed pLGIC hydrophobic gate requires too high a voltage to occur physiologically but may inform designs for switchable nanopores. Global analysis of ∼200 channels yields a simple heuristic for structure-based prediction of (closed) hydrophobic gates. Simulation-based analysis is shown to provide an aid to interpretation of functional states of new channel structures. These studies indicate the importance of understanding the behavior of water and ions within the nanoconfined environment presented by ion channels.

Introduction: Ion Channels and Hydrophobic Gating

Ion channels are proteins that form nanoscale pores in biological membranes. Channels play a central role in many aspects of cell biology and physiology, especially in excitable cells of the nervous system, and consequently are a major class of drug targets. Advances in structural biology, most recently in cryo-electron microscopy (cryo-EM), have revealed over 1000 three-dimensional structures of ∼160 different ion channel proteins. This highlights the challenges of relating the molecular structure of channels to their physiological function. A key aspect of ion channel function is gating, namely, the controlled switching of a channel between a closed, i.e., non-conductive, state and an open, i.e., conductive, state. Gating may be controlled by a number of physiological factors, including the voltage difference across a cell membrane, membrane stretch, and the binding of specific ligands to receptor domains of a channel protein.

At the level of the transmembrane nanopore formed by a channel protein, a number of possible mechanisms exist for switching between a closed and an open state. One mechanism, hydrophobic gating, is of particular interest from the perspective of physical chemistry, as it exploits the unique properties of water nanoconfined in a hydrophobic pore. When in a constricted (radius <0.5 nm) hydrophobic pore, water may exist in a vapor state, such that the pore is dewetted, thus preventing the permeation of ions (Figure A). Such a pore is therefore functionally closed, even when it is not physically occluded. This hydrophobic gating mechanism has been evaluated using simple computational models (Figure B). More recently, a number of channel structures have suggested that this mechanism may operate in several ion channel species. In particular, the pentameric ligand gated ion channels (pLGICs[1])—neurotransmitter receptors which play key roles in the nervous system—appear to be controlled via hydrophobic gating (Figure C).[2] Alongside other channels with a hydrophobic gate (e.g., MscS[3−5] and BEST[6,7]), the pLGICs provide a test bed for studies of hydrophobic gating in ion channels more generally. In this review, we will focus on the use of molecular dynamics (MD) simulations to complement experimental structural advances by providing a powerful approach for dissecting the mechanism of hydrophobic gating and aiding functional assignment of novel channel structures.

Figure 1

Hydrophobic gating in ion channels and selected simulation systems. (A) Schematic of hydrophobic gating, showing the transition between a closed state, dewetted (white) in a hydrophobic constriction in the center of the channel (dark gray) and an open state in which the channel is hydrated throughout. The lipid bilayer is shown in pale gray, water in cyan, and ions in blue and red. Figure adapted with permission from ref (74). Copyright 2019 National Academy of Sciences. (B) Simple model of a nanopore (cyan), embedded in a membrane-mimetic slab (gold), with water molecules (red/white) on either side and within the pore. Figure reproduced with permission from ref (18). Copyright 2003 Wiley. (C) Structure of a pentameric ligand gated ion channel (pLGIC), the 5-HT3 receptor in an open state, shown as a Cα trace with the pore-lining surface (as determined by CHAP) colored on hydrophobicity. The arrow indicates the location of the hydrophobic constriction/gate formed by the ring of L9′ residues. Figure adapted with permission from ref (125). Copyright 2019 Elsevier. (D, E) Complete channel protein (D) compared with a model nanopore (E) corresponding to the M25 pore-lining helix bundle, both illustrated for the 5-HT3 receptor, with the lipid bilayer in gray/brown. Figures reproduced with permission from ref (68). Copyright 2020 American Chemical Society.

Computational Methods and Models

MD simulations have been widely employed to study hydrophobic gating, both in simple model systems (Figure B and below) and in biological ion channels. Studies of large and complex ion channel structures (Figure C) either may employ the intact channel protein embedded in a lipid bilayer (Figure D) or may focus on either the transmembrane (TM) domain or the immediate pore-lining structure (Figure D,E). In the case of the pLGICs, the pore-lining structure consists of a bundle of five M2 α-helices (M25). This provides a convenient model of a biological nanopore of sufficient complexity to merit in depth examination, while sufficiently simple (and small) to enable the application of advanced and hence more computationally demanding methodologies.

MD simulations may be used to study the dynamic interactions of water and of ions with different conformational states of an ion channel, both closed and open. Potential of mean force (PMF) calculations provide estimates of the free energy landscape encountered by ions as they pass through a transmembrane pore. Computational electrophysiology (Comp ePhys) methods,[8−11] in which a transmembrane potential is imposed upon a membrane-embedded channel or pore, may be used to simulate permeation of ions, thereby enabling direct comparison with experimental single-channel conductance measurements. There have been a number of authoritative general reviews of simulation approaches to ion channels.[12−16] Here we will focus on simulation studies of hydrophobic gating and their relationship to conductance properties of some recently determined channel structures.

Simulations and Hydrophobic Gating

Early studies focused on the behavior of water confined within model nanopores (Figure )[17−26] and in carbon nanotubes (CNTs),[27] thereby establishing the general features of hydrophobic gating. These studies demonstrated that water confined within a narrow (radius <0.5 nm) hydrophobic pore can coexist in a vapor and liquid state, with the free energy difference between these states determined by the pore radius and the polarity of the pore lining (Figure ) and also by the voltage across the pore/membrane system. Thus, a hydrophobic nanopore could be switched from a dewetted (i.e., closed) state to a hydrated (i.e., open) state either by a small increase in the pore radius and/or the pore polarity or by imposition of a relatively high (i.e., greater than physiological) voltage difference across the pore. These theoretical considerations informed, e.g., the design of synthetic nanopores which exhibited voltage-sensitive hydrophobic gating.[28] However, a paucity of high-resolution structures precluded wider evaluation of the significance of hydrophobic gating for biological ion channels.

Figure 2

Hydrophobic gating in a simple model nanopore. (A) Pore model consisting of methane-like particles (dark gray) of van der Waals radius 0.195 nm arranged to form a membrane-like slab containing a hydrophobic pore of minimum radius R. A water molecule (w) is shown drawn to scale. Figure adapted with permission from ref (19). Copyright 2003 National Academy of Sciences. (B) Hydration probability (ω) vs pore radius for model nanopores. The gray region indicates radii smaller than the radius of a water molecule (0.14 nm). Data points are obtained from MD simulations with the errors estimated from block averages. The continuous lines are fits of a simple model to data points for hydrophobic (black), amphipathic (red), and polar (green) pores. The vertical line indicates the approximate radius of the closed state of a pLGIC. (C) Water and ions in a hydrophobic nanopore. Liquid–vapor oscillations of water shown as the hydration probability ω in a R = 0.65 nm hydrophobic nanopore (top panel) in the presence of a 1.3 M NaCl solution. As indicated by the number N of ions within the pore, Na+ (middle) and Cl– ions (bottom) are only observed in the pore when there is also liquid water (ω ≈ 0.8) present. Permeation events are indicated by triangles; ions do not permeate the pore during the vapor phases. Parts B and C were reproduced with permission from ref (20). Copyright 2004 IOP Publishing.

Hydrophobic Gating in Ion Channels

Several simulation studies explored hydrophobic gating using early structures of pLGICs, including a low-resolution structure of a nicotinic acetylcholine receptor (nAChR[29,30]), and an X-ray structure of GLIC, a bacterial pLGIC.[31,32] Hydrophobic gating was also suggested by simulations of the bacterial mechanosensitive channel MscS.[3] In each case, a closed state of the channel was identified in which the hydrophobic gate region of the TM domain was dewetted.

Hydrophobic gating, and more generally speaking hydrophobic constrictions contributing to gating, have since been invoked for a wide range of channels (recently reviewed in ref (33)), including, e.g., potassium channels (Kv channels,[34] K2P channels,[35] BK channels,[36] and NaK2K channels[37]), Hv (voltage gated proton permeable) channels,[38,39] the CorA magnesium channel,[40] and Orai.[41] Hydrophobic gate-like structures have also been suggested for various transporters including, e.g., the ABCG2 multidrug transporter[42] and the NhaP Na+/H+ antiporter.[43] Extension of the hydrophobic gating concept to transporters is of interest given the proposed role of water filled channel-like states in vSGLT and related transporters.[44] Over the past few years, multiple structures of bacterial and animal pLGICs have been determined, for which the corresponding functional state (e.g., open vs closed) is not always clear-cut. In this situation, MD simulations may be used to aid the annotation of an ion channel structure. An example of this is provided by a crystal structure (PDB id 4PIR)[45] of the 5-HT3 receptor (5-HT3R), a biomedically important member of the pLGIC family present within the mammalian nervous system which is activated by the neurotransmitter 5-hydroxytryptamine (5-HT aka serotonin). At the time of the structure determination, the authors stated “the 5-HT3 receptor pore state is not clearly defined”.[45] Simulations of the corresponding M25 pore-lining domain (Figure A) in a phosphatidylcholine (PC) bilayer, using the TIP4P water model, revealed local dewetting (Figure B) of the pore in the vicinity of a ring of hydrophobic leucine (L9′) side chains which form the hydrophobic gate.[46] Free energy landscape (i.e., PMF) calculations revealed an energetic barrier in this region, both for water and for ions (Figure C). Extended (microsecond) MD simulations of the intact 5-HT3R in a phospholipid bilayer also demonstrated local dewetting of the hydrophobic gate, further supporting the assignment of the structure to a closed state of the channel.[47]

Figure 3

Analysis of hydrophobic gating for a closed state of a pLGIC (the 5-HT3R, PDB id 4PIR). Figures modified with permission from ref (46). Copyright 2016 Cell Press. (A) M2 helices lining the pore, showing two of the five pore-lining helices along with the pore surface (gray). Pore-lining side chains are shown, with the hydrophobic gate formed by the L9′ and V13′ rings in blue. (B) Positions of water molecules projected onto the pore (z axis) as a function of time for an MD simulation of the M25 helix bundle shown in part A embedded in a PC bilayer. (The dashed horizontal lines indicate the phosphate headgroup positions of the lipid bilayer.) Each water molecule is represented by a light blue circle. The trajectories of a number of individual water molecules are illustrated using darker blue. Thus, the intermittent white region around z = 0 corresponds to the dewetted L9′ hydrophobic gate. (C) Free energy profiles for water and ions along the axis of the 5-HT3R (PDB id 4PIR) pore. These were estimated as potentials of mean force (PMF) for single ions (Na+, red; Cl–, green) or single water molecules (blue) as a function of position along the pore (z) axis. The gray shading represents the extent of the pore, with vertical lines indicating the positions of the pore-lining side chains.

Computational Physical Chemistry of Hydrophobic Gating in pLGIC Ion Channels

The pLGICs have been used to explore the computational physical chemistry of water in ion channels. The results from these studies of complex biological channels may be compared with a large body of studies of water in nanopores, the latter ranging from simplified (conceptual) models to CNT porins and synthetic biomimetic channels.[48] In particular, it is important that advances in modeling water in nanopores[48] are applied to simulations of hydrophobic gating, especially models of molecular polarizability.[49−52]

M25 nanopore models of pLGICs have been used to explore the sensitivity of simulations of hydrophobic gating to the water model employed. Initial simulations both of channel structures and of simple models used “tried and tested” water models (e.g., TIP3P or SPC). Comparative studies[53] have suggested that some more recent water models (e.g., TIP4P/2005[54] and OPC[55,56]) exhibit better quantitative agreement with experimental data for, e.g., interfacial properties of water and/or interactions with ions. Simulations of M25 models derived from early cryo-EM structures of the glycine receptor (GlyR),[57] a biomedically important anion selective pLGIC, suggested that differences between water models with respect to hydrophobic gating were relatively small, at least with respect to a closed/desensitized state of this channel (PDB id 3JAF) when compared for the TIP3P, SPC/E, and TIP4P models.[46] This comparison has recently been extended, using M25 nanopore systems from four conformations of the 5-HT3 receptor, and comparing both additive (i.e., fixed charge) and polarizable models of water.[58] For a closed state (PDB id 4PIR; see above) of the 5-HT3R, again only minor differences in the behavior of water were seen (Figure ). In contrast, for an open-state structure (PDB id 6DG8), there were clear dependencies of the wetting/dewetting profile along the length of the pore on the water model used (Figure B). Thus, the pore was incompletely hydrated (i.e., partially dewetted) when using older additive models for water (e.g., TIP3P and SPC/E) while more fully wetted with either recent additive models (e.g., TIP4P/2005[54] and OPC[55,56]) or with the AMOEBA[59] polarizable model. This suggests that for a hydrophobic gate structure which is on the “edge” of wetting the choice of water model in simulations may be critical.

Figure 4

Simulations of wetting/dewetting and sensitivity to water models. (A) Dewetted pore from a simulation of the closed state 5-HT3R pore (PDB id 4PIR) using the mTIP3P water model. Protein and lipid molecules are in gray and waters in blue. The red circle indicates the dewetted region of the pore. (B) Water models compared for a closed state (PDB id 4PIR) and an open-state (PDB id 6DG8) structure of the 5-HT3R pore. Water density profiles along the pore axis are shown for the mTIP3P (green), TIP4P/2005 (yellow), AMOEBA03 (red), and AMOEBA14 (black) water models. The dashed horizontal line indicates the density of bulk water, and the dashed vertical lines and the shaded background denote the extent of the protein and the hydrophobic gate region, respectively. Single-letter codes for the pore-lining amino acid side chains are given at the top of each panel. Figure modified with permission from ref (58). Copyright 2020 American Chemical Society.

The use of a polarizable force field (as opposed to fixed point charge additive models) also influences the energetics of interactions of ions and water with the pore, as can be seen from free energy profiles (as evaluated by PMFs) for a single ion (e.g., Na+ or Cl–) moved through a pore embedded in a lipid bilayer. This is illustrated for the open-state (6DG8) M25 pore of the 5-HT3 receptor in Figure A.[58] Comparing single-ion PMFs estimated using an additive force field (CHARMM36m with TIP3P water) with those obtained using a polarizable force field (AMOEBA14), it can be seen that the overall shape and features of the profiles are preserved. Thus, for both Na+ and Cl– ions, there is an overall barrier to be crossed. This barrier is lower for Na+ ions (the 5-HT3R is cation selective) than it is for Cl– ions. However, the polarizable force field results in more rugged energy landscapes for either species of ion, suggesting that there may be preferential interactions in regions along the pore where the ion can induce dipoles in neighboring protein atoms.

Figure 5

Polarizability and interactions with ions. (A) Single-ion free energy (i.e., PMF) profiles for Cl– (red) or Na+ (black) ions in the pore of the open-state 5-HT3R (PDB id 6DG8). The vertical dashed lines denote the extent of the protein, and the yellow bar represents the position of the lipid bilayer. The hydrophobic gate region is denoted by the gray background shading. Single-letter codes for the pore-lining amino acid side chains are given at the top of each panel. The collective variable (CV) is defined as the distance along the z-axis between the ion and the protein center of mass and thus is zero at the center of the pore. Profiles are in each case shown for an additive (CHARMM36m/mTIP3P) and for a polarizable (AMOEBA14) force field. (B) Cl– interactions with the hydrophobic surface lining the pore when the polarizable AMOEBA14 force field is used. A snapshot from an umbrella sampling window is shown. The Cl– ion is shown as a yellow van der Waals sphere, while oxygens of water molecules in the first and second hydration shells are shown in cyan and blue, respectively. A zoomed-in image (right) of a Cl– ion interacting with three hydrophobic side chains (L9′, V13′, I17′) is shown. Figure modified with permission from ref (58). Copyright 2020 American Chemical Society.

Interestingly, polarizable Na+ and Cl– ions exhibit rather different behaviors close to the hydrophobic pore-lining side chains in the gate region of the 5-HT3R pore. Thus, a Cl– ion relinquishes part of its hydration shell to form a close association with the hydrophobic surface, whereas a Na+ ion retains most of its hydration shell. This difference between Na+ and Cl– ions at hydrophobic pore surfaces is not seen for additive force fields. Interaction of Cl– ions with hydrophobic residues can be seen at a number of locations along the pore, such that at one position a Cl– ion is hydrated by 4–5 inner-shell waters on one side of the ion but on its other side ∼2 waters have been displaced and instead the ion interacts directly with a hydrophobic pore surface (Figure B). This suggests that, when a polarizable force field is used, Cl– ions may form favorable interactions close to a hydrophobic surface of the pore, which agrees with, e.g., simulation studies of halide ions at water/vapor[60−62] and at water/hydrophobic interfaces in general.[63] Interactions with hydrophobic (alkyl) groups have been observed in structures of chloride selective anionophores (e.g., biotin[6]uril) esters) and of chloride ion transporter protein structures (e.g., the chloride-pumping rhodopsin from Nonlabens marinus(64)). These subtle differences between additive and polarizable force fields are therefore likely to prove important when estimating the energetics of ion permeation through pLGICs (see, e.g., refs (65 and 66) for recent applications).

An important but potentially underexplored aspect of hydrophobic gating is the effect of a transmembrane electric field in promoting wetting of an otherwise dewetted hydrophobic gate. Initially explored in simple model systems (e.g., refs (24 and 25) also recently reviewed in ref (48)), this has also been observed in the bacterial mechanosensitive channel MscS[4] and in model protein nanopores.[67] In these systems, the presence of a strong (e.g., greater than physiological) electric field corresponding to a transmembrane voltage difference of ΔV = 0.5–1 V drove hydration of an otherwise dewetted hydrophobic gate. The 5-HT3R M25 system has been employed as a biologically realistic model nanopore[68] which enables detailed exploration of this behavior, including the sensitivity to the water model employed (Figure ). It can be seen that complete wetting/hydration of the closed state (PDB id 4PIR; see above and Figures and 4) requires a field strength in excess of 100 mV/nm (corresponding in this system to a transmembrane voltage difference of ΔV = 0.85 V) when the TIP3P water model is used. This threshold depends on the water model and is increased to ΔV > 1 V for TIP4P/2005. Electric field dependent wetting of a nanopore is well described by a simple model embodied in the following equationwhere ⟨ω⟩ = the time-averaged hydration probability, ΔΩ0 = ΩV – ΩL (i.e., the difference between the free energies of the liquid and vapor states in the absence of an E-field), and β = 1/kBT. When ΔΩ0 < 0, the pore is hydrophobic and favors a vapor state. The effect of the electric field, E, is represented by a second free energy term, in which m denotes the strength of the coupling between the hydration probability and the magnitude of the field and where EINT accounts for the horizontal offset of the response curve by an intrinsic electric field arising from the nanopore structure. Thus, m represents the difference in ability to store electrical energy between a water-filled (high-dielectric) space vs an empty (low-dielectric) space. This in turn is related to the wettable volume of the pore and to the local dielectric constant of the water model. As can be seen from Figure D, this model fits simulation data for the effect of the applied electric field on hydration probability for both closed (PDB id 4PIR) and open (PDB id 6DG8) states of the 5-HT3R M25 model pore. Although, given the high fields involved, such effects are unlikely to result in a switch between dry/wet and closed/open states of biological ion channels under physiological conditions, they may provide the basis of electric field switchable biomimetic pores, as has been demonstrated experimentally.[28,69]

Figure 6

Electrowetting of the hydrophobic gate in a nanopore. (A) Simulation system: an M2 helix bundle nanopore (green) is embedded in a lipid bilayer (brown). Water is represented as a transparent surface, and Na+ and Cl– ions are shown in red and yellow, respectively. The effect of a transmembrane potential is modeled by applying a constant electric field to atomic charges in the system. The transmembrane voltage is reported as VIC– VEC, where IC = intracellular (negative z) and EC = extracellular (positive z). (B) Time series of hydration probability (ω) at three different electric field strengths. The gray line in the background represents the water density in the hydrophobic gate normalized to the density of bulk water. The colored lines represent discretization of this via a threshold crossing algorithm (see part D below). (C) Time-averaged water density profiles as a function of electric field strength (E). The shaded background and vertical dashed lines indicate the extent of the hydrophobic gate and of the transmembrane domain, respectively. The horizontal dashed line represents the density of bulk water. (D) Fitting a simple model to the simulation data for hydration probability ⟨ω⟩ as a function of the external electric field, E. Data points and error bars represent the mean hydration probability and its standard error over three independent simulations. The solid lines are from fitting a nonlinear model (see equation and main text). Simulations were performed on the M2 helix nanopore from the closed (PDB id 4PIR) and open (PDB id 6DG8) states of the 5-HT3 receptor. The data shown in parts B and C are based on simulations of the M2 helix nanopore using the mTIP3P water model. Figure modified with permission from ref (68). Copyright 2020 American Chemical Society.

Global Approaches to Prediction of Hydrophobic Gating

MD simulations have also been used to explore hydrophobic gating in a range of ion channel proteins other than pLGICs, including bestrophin (BEST1[70,71]) and TMEM175,[72] both of which have a putative gate formed by three adjacent rings of hydrophobic residues.[6,7] A hydrophobic gate has also recently been suggested by simulations of the two pore channel TPCs.[73] In BEST1, the three rings of hydrophobic residues forming the gate are I76, F80, and F84. These have been mutated in silico to explore the relationship between the nature of the amino acid side chains forming the gate and the height of the energetic barrier for water crossing that gate. Replacement of the I76, F80, and F84 rings by three rings of identical hydrophobic aliphatic side chains (e.g., III, LLL, or VVV) maintains the barrier to water (and hence ion) permeation, whereas replacement with a polar side chain (e.g., T, threonine) enables the pore to wet. Thus, the pore radius profile at the gate is approximately identical for VVV and TTT (the side chains of V and T are isosteric), but the gate is dewetted in the former case while fully hydrated in the latter case. This result suggests the need for a more wide-ranging examination of hydrophobic gates in ion channels of known structures.

A global approach to hydrophobic gating in ion channels analyzed simulations of 190 different channel structures in order to define the relationship between the local radius and hydrophobicity of the channel and the height of the resultant free energy barrier to water permeation at the gate (Figure A,B).[74] A clear-cut relationship was demonstrated, which in turn enabled development of a simple heuristic for predicting whether a channel is likely to be closed or open based simply on the radius and hydrophobicity of the pore/gate region. The application of this heuristic for, e.g., two recently determined structures of a closed (PDB id 6V4S) and an open (PDB id 6V4A) state of DeCLIC, a bacterial pLGIC,[75] is shown in Figure C. This structure-based method has been shown statistically to perform better than, e.g., a prediction based on pore radius profile alone.[74] Importantly, this analysis demonstrates that to a first approximation the behavior of a hydrophobic gate can be described well by just the local pore radius profile and hydrophobicity, which enables robust prediction of the functional state of new channel structures and also provides a clear design principle for hydrophobic gates in synthetic nanopores.[76]

Figure 7

A global survey of hydrophobic gating in ion channels. (A) Montage of the structures of 190 channels surveyed, each viewed down the pore axis. (B) Schematic of hydrophobic gating as a function of (hydrophobicity, radius) of the transmembrane pore. The surface shows the free energy of water within a channel as a function of (hydrophobicity, radius) corresponding to the data set of unique channel structures in part A. Schematic depictions of dewetted (closed) and hydrated (open) states of channels are shown for the two main regions of the data. (C) Illustration of the heuristic method for two structures of the bacterial pLGIC DeCLIC, in the presence and in the absence of Ca2+, corresponding to closed (PDB id 6V4S) and open (PDB id 6V4A) states, respectively. The pore-lining surfaces (colored on hydrophobicity, pale-brown corresponding to maximum hydrophobicity) are shown. In the graphs, for each pore-lining side chain, the channel pore radius at the residue is plotted against the corresponding local hydrophobicity value. The sum of shortest distances between the dashed (1 RT) contour line and all points falling below it (colored red) are used as a score for identifying closed gates. A structure is predicted to be in a non-conductive state if it has a value of ∑d > 0.55. Figures modified with permission from ref (74). Copyright 2019 National Academy of Sciences.

Annotating New Structures: Three Recent Case Studies

A simulation-based approach has been used for the functional annotation of new structures of pLGICs and of other ion channels. Two recent studies have used simulations of the behavior of water within the TM pore domain to assign functional states to newly determined cryo-EM structures for the 5-HT3R corresponding to the channel in closed vs open conformations. These simulations and structures[65,77] added to the two previous studies of the closed state (PDB id 4PIR) as discussed above.[46,47] Two independent but parallel studies identified an open state (PDB ids 6DG8 and 6HIN, respectively) which was fully permeable to water and ions and also “preopen” (or perhaps desensitized) states (PDB ids 6DG7 and 6HIO, respectively) which remained dewetted in the hydrophobic pore region. Thus, simulations of water behavior have contributed to our understanding of the functional importance of the different conformational states determined for the 5-HT3R. A comparable situation has recently been seen for a bacterial mechanosensitive channel, YnaI, for which there is a closed-like low-conductance state (PDB id 6ZYD) which in simulations exhibits a partially dewetted hydrophobic gate (with a free energy barrier of ∼1.5 kJ/mol, i.e., 0.6 RT to water permeation) and also a fully open hydrated state (PDB id 6ZYE) which presents no barrier to water permeation.[78]

A more complicated case of assignment of functional states to structures is presented by the glycine receptor (GlyR; discussed briefly above). A landmark cryo-EM study revealed structures of detergent-solubilized zebrafish GlyR, in both closed and open states.[57] The nature of the hydrophobic gate in these three structures (PDB ids 3JAD, 3JAE, and 3JAF for the closed, open, and desensitized states, respectively) was analyzed by simulations, and as noted above, the desensitized state (3JAF) was used for an initial exploration of the robustness of hydrophobic gating simulations to the water model employed.[46] More recently, there has been a discussion of the nature of the open-state structure determined using detergent-solubilized GlyR protein, focusing on whether or not this may represent as a “superopen” state, the relationship of which to the physiological open state in a cell membrane (i.e., a lipid bilayer) is unclear.[79−81] There has been an attempt to resolve this via the use of Comp ePhys simulations to predict ionic conductances (see next section) for comparison with experimental single-channel measurements. However, the results of this comparison and their interpretation have proved to be somewhat controversial and are the subject of an ongoing discussion.[82,83]

More recently, cryo-EM has been used to determine multiple structures of the GlyR in lipid nanodiscs, which are believed to provide an environment more closely resembling the lipid bilayer in a native cell membrane.[84,85] In both studies, MD simulations were used to aid identification of a physiologically open state. Here we will focus on the use of three different levels of computational study to explore the nature of the hydrophobic gate in these bilayer-embedded GlyR structures.[85] In particular, we will compare the closed (apo) state of the GlyR (PDB id 6UBS) with an open state of the channel (PDB id 6UD3) obtained by studying the GlyR in a complex with both the agonist (i.e., activator) glycine plus an open channel blocker, picrotoxin (PTX).

Comparison of the pore-lining surfaces of the 6UBS (apo) and 6UD3 (+Gly, +PTX) states of the GlyR (Figure A) reveals differences in the radius profiles of the transmembrane pore, especially in the region of the hydrophobic gate formed by the key ring of L9′ side chains. This is especially evident if the structure-based heuristic procedure (see above) is adopted, displaying the positions of the pore-lining side chains in the local (hydrophobicity, radius) plane (Figure B). For the 6UBS (apo) state, there are multiple residues below the heuristic cutoff line, providing an unambiguous prediction that this is a closed state of the channel. In contrast, there are no points below the line for the 6UD3 structure, predicting this to be a fully open state. This is further confirmed by short simulations (Figure C) which demonstrate an energetic barrier to water (and hence by proxy to ions) in the region of the L9′ gate for the 6UBS structure, whereas for the 6UD3 structure there is no barrier for water permeation along the length of the transmembrane pore. Finally, Comp ePhys simulations on the intact receptor molecule embedded in a phospholipid bilayer in the presence of an externally applied electrostatic field (corresponding to 500 mV) revealed that, while Cl– ions passed through the transmembrane pore in the 6UD3 conformation, no ions passed through the pore for the 6UBS conformation, the L9′ gate region of which remained dewetted even in the presence of a transmembrane voltage (Figure D). Thus, we can see how structural/heuristic analysis derived ultimately from a global simulation study of hydrophobic gates in ion channels can be used to bring clarity to a complex situation with multiple states of a pLGIC and how this initial analysis may then be validated and extended by further simulations of the behavior of water and ions in the vicinity of a hydrophobic gate.

Figure 8

Application of three levels of simulation analysis to novel pLGIC channel structures. Figures modified with permission from ref (85). Copyright 2020 Nature Research. (A) Two structures of the GlyR, in a closed (apo; PDB id 6UBS) and an open (GlyR + glycine + PTX; PDB id 6UD3) state. For each structure, the pore-lining surface is colored by hydrophobicity (green for hydrophilic to brown for hydrophobic) as estimated by CHAP.[125] (B) Likelihood of pore closure by an energetic barrier corresponding to dewetting at a hydrophobic constriction, evaluated according to a heuristic method based on simulation of water behavior in 190 ion channel structures (see Figure ).[74] Pore-lining side chains are indicated using their local pore hydrophobicity and radius as coordinates. The subset of points falling below the dashed classification line is used to calculate a heuristic score. A cutoff of ∑d > 0.55 predicts that a channel structure contains a hydrophobic barrier to water and ion permeation. (C) Water free energy profiles for the two GlyR structures. The free energy profiles in dark blue were derived from three independent 30 ns simulations during which positional restraints were applied to protein backbone atoms. Time-averaged profiles derived in each case from the final 20 ns of 200 ns unrestrained simulations are shown in gray. (D) Water (pale blue) and ion (Cl–, mid-blue; Na+, dark blue) trajectories projected onto the pore (z) axis for the two GlyR structures. The L9′ residues are located at z = 0 nm. Simulations of the GlyR in a PC bilayer were in 0.5 M NaCl with a transmembrane potential of +500 mV applied via a uniform external electric field (with positive potential on the cytoplasmic, i.e., negative z side). Positional restraints were applied to protein backbone atoms, in order to preserve the experimental conformational state while permitting rotameric flexibility of amino acid side chains.

Predicting Conductance

From a physical chemistry perspective, ion channels provide an ideal opportunity to relate simulation studies of water and ion movement directly to single-molecule (i.e., patch clamp) measurements of channel conductance. This has been explored in some detail for potassium channels, the exquisite selectivity of which reflects ion permeation in a partly or completely dehydrated state through a conformationally flexible filter. However, there is an ongoing debate over the exact mechanism of potassium channel permeation,[86−92] which makes quantitative comparison of simulated ion permeation rates and experimental single-channel conductance challenging. For pLGICs, which when open are thought to allow selected ions to permeate with their first hydration shell intact, one might anticipate that it should be easier to match experiments and simulation.[79] However, recent studies highlight both computational and experimental challenges which remain to be addressed before we can accurately predict the rate of movement of hydrated ion flow through a nanopore of known structure and dimensions.

To explore in more detail the relationship between pore hydrophobicity, geometry, and channel conductance, we selected a subset of 22 channel structures. These were of channel proteins that are conductive of hydrated ions, possess hydrophobic gates, and have more than one structural state determined at a reasonable resolution. These channel structures (pLGICs: 5HT3R, GABAAR, GlyR, nAChR, GLIC; and also MscL and Orai) were selected such that the pore radius at their water free energy maximum is above 0.13 nm; i.e., they are not physically occluded. Comp ePhys simulations were run for each of these channel structures, in the presence of 0.5 M NaCl, and at ΔV = 500 mV. Backbone restraints were applied, so as to prevent conformational drift from the experimentally determined structures, and three repeats of 200 ns were performed in each case. Shorter (3 × 30 ns) equilibrium (i.e., no ΔV imposed) simulations were run in order to obtain water densities along the length of the transmembrane pore. Further details are provided in Supporting Information Table S1.

The results are summarized in Figure , which shows simulated conductance vs minimum water density within the transmembrane pore for these 22 channel structures. Below a minimum water density of 12 nm–3 (corresponding to 36% of bulk water density) within a channel, it appears to be functionally closed, i.e., non-conductive. Above this water density, either high conductance (i.e., open) or low but finite conductance (often corresponding to a desensitized or intermediate state) channels are generally seen. The 12 nm–3 water density cutoff seems to be a reliable predictor of a closed state of a pLGIC. However, further data (i.e., structures) are needed to obtain a clear distinction between an open (high conductance) and a desensitized state, as both may exhibit a water density >12 nm–3 within the pore. The key structural difference may lie in the disposition of charged pore-lining side chains at the non-hydrophobic mouth of the pore, but it is difficult to be certain on the basis of the current available structures.

Figure 9

Pore hydration level as an indicator of ion channel conduction state. Data from analysis of Comp ePhys simulations for 22 channel structures (pLGICs: 5HT3R, GABAAR, GlyR, nAChR, GLIC; with also MscL and Orai), selected such that the pore radius at their free energy maxima is above 0.13 nm, i.e., they are not physically occluded. Simulated ion conductances vs minimum water densities for the 22 selected channel structures. Conductance and density values are, respectively, averaged between triplicates of 200 ns Comp ePhys simulations (in the presence of 0.5 M NaCl) and 30 ns equilibrium simulations. Red, gray-green, and blue points correspond to closed, desensitized/intermediate, and open-state structures, respectively. The gray vertical dotted line represents a bulk water density of 33.4 nm–3. The red dotted line represents a heuristic cutoff of 12 nm–3.

Future Challenges and Outlook

The studies described above have largely focused on the interactions with water and ions in relation to the conductance properties of an ion channel, e.g., closed/non-conductive vs open/conductive. An ongoing challenge is to understand in detail the pathways between these different conformational states and the mechanisms controlling the transitions between them. For the pLGICs, this centers around ligand induced conformational changes of the extracellular domain and how these are transmitted to changes in the conformational state of the transmembrane pore. This has been the subject of a number of simulation studies (see, e.g., refs (93 and 94)) and remains the focus of several ongoing studies, and so we will only address this briefly here. The main challenge is the time scale of conformational transitions associated with channel gating (i.e., msec[95,96]), which is currently difficult to address directly via MD simulations. A range of approaches have been applied including direct (long) simulations of 5-HT3R[65,97] and other pLGICs, e.g., GluCl,[98] and the use of a range of enhanced sampling approaches,[99] including, e.g., the use of strings and swarms of simulations[100,101] and of extensive cloud-based nonequilibrium simulations.[102] A promising approach is to combine extended MD simulations initiated from multiple starting states as determined by, e.g., cryo-EM (seen, e.g, in a recent study of the bacterial channel GLIC[103]) or NMR (as has recently been applied to gating modes of potassium channels[92]). It is likely that such an approach, combined with, e.g., Markov state modeling,[104−107] will enable a rigorous integration of insights from structure-based simulations and from single-channel kinetic analysis and modeling.

In addition to computational approaches, advances in our understanding will also be driven by higher resolution and more complete structural data for pLGICs and related ion channels. Recent improvements in the resolution (to 1.7 Å) of membrane (ion channel) protein cryo-EM have revealed individual water molecules in pLGIC structures.[108] Elements other than proteins, e.g., covalently bound glycans and non-covalently bound lipids,[109] can also be revealed by cryo-EM, and in both cases, simulations will be key to understanding the role of these interactions in the functional properties of pLGICs.[110,111] In particular, structural and simulation studies together can provide insights into the functional roles of specific lipid molecules bound to ion channels.[112]

It would of course be highly desirable to experimentally validate hydrophobic gating of ion channels. This could be approached using, e.g., time-resolved FTIR spectroscopy, which has been used to determine the configuration of water molecules in bacteriorhodopsin[113] and in halorhodopsin.[114]

Alongside advances in structural data, improved computational models continue to be developed. In particular, more widespread use of polarizable force fields[51] and their efficient computational implementation[115] promise to refine our models of channel/water/ion interactions.[51] It is therefore important to consider whether we observe meaningful improvements in our understanding of channel function given the increased computational cost of using polarizable force fields. The increased “ruggedness” of free energy landscapes for ion permeations obtained using polarizable force fields (see above; Figure ) is likely to influence ion movement within a pore and thus to change predictive estimates of ion conductance. Furthermore, the behavior of anions at water/hydrophobic interfaces is known to be sensitive to inclusion of polarizability.[60,63] Polarizability is also likely to be of considerable importance for accurate modeling of divalent cations, e.g., Ca2+ ions and their interactions.[116] An area which merits further exploration in terms of applications to ion channels is the use of charge scaling, i.e., electronic polarizability treated implicitly via the electronic continuum correction (ECC) model.[117,118] This has been shown to accurately mimic polarizability in terms of preferential localization of anions at water/hydrophobic interfaces[63] (as is also seen in ab initio MD[119]) and has been used for, e.g., Ca2+ interactions with anionic lipid bilayers.[120] This approach therefore offers the possibility of increased accuracy in ion channel simulations at little extra computational cost. However, as yet, the ECC model has not been tested for (membrane) proteins.[52]

In combination with implementation of biomolecular simulation codes on exascale computing resources[121−124] and growth in the number of ion channel structures,[125] these advances in force fields for channel simulations offer the possibility of, e.g., global comparisons of ion channel conductance behaviors obtained by large scale Comp ePhys surveys with more accurate models of ion/water/channel interactions, thus yielding new insights into the atomic level relationship between structure and physiological function in biological ion channels.