Identification of Hydrophobic Gates by Simulation-Enabled Analysis of the Ion Channel Structural Proteome

Ion channel proteins control ionic flux across biological membranes through conformational changes in their transmembrane pores. An increasing number of channel structures captured in different conformational states are being determined as a result of advances in membrane protein structural biology. However, newly resolved structures are commonly classified as either open or closed based solely on the physical dimensions of their pore and it is now known that more accurate annotation of their conductive state requires an additional assessment of the effect of pore hydrophobicity. A narrow hydrophobic gate region may disfavour liquid-phase water, leading to local de-wetting and an energetic barrier to water and ion permeation without steric occlusion of the pore. Here we quantify the combined influence of radius and hydrophobicity on pore de-wetting by applying molecular dynamics-based calculations to nearly 200 ion channel structures. Using machine learning, we propose a simple heuristic model that allows the rapid prediction of the state of hydrophobic gates based on analysis of both local radius and hydrophobicity in ion channel structures and nanopores.


Introduction
Ion channel proteins are water-filled, ion-conducting pores that are key components of biological membranes (Hille, 2001). In a resting or inactivated state, the flow of ions through the pore may be interrupted at one or more positions along the channel, termed gates. Upon activation, for example in response to ligand binding or a change in the membrane potential, conformational changes generally lead to expansion of the channel pore at its gate region(s) thereby switching it from a closed (i.e. non-conducting) to open (conductive) conformation. Ion channels represent attractive therapeutic targets and so there is considerable interest in elucidating the mechanisms underlying this process in many different ion channel families. In most cases, the determination of channel structures in various different conformational states provides the key molecular basis for understanding how this 'gating' is regulated, often supported by insights from electrophysiological and other biophysical approaches.
For a newly determined ion channel structure, its conductive state is most frequently inferred by measuring the physical dimensions of its transmembrane pore, displayed as a profile of pore radius along the central axis along with an image of the pore-lining surface. There are several methods for this, with one of the most widely used being the HOLE program (Smart et al., 1996). Such pore radius profiles provide an indication of the maximum size of ions that might be accommodated within a transmembrane pore and steric constrictions narrower than the radius of a water molecule (~0.15 nm) are considered as potential gates.
However, the permeation of ions and water through a given region in a sub-nanometre pore is influenced not only by its radius, but also by the local hydrophobicity of the pore lining. Ion permeation may readily occur through polar regions only just larger than the hydrated radius of the permeating species, but a hydrophobic pore segment of comparable dimensions may undergo spontaneous de-wetting to form a local nanoscale region devoid of water and ions (Beckstein et al., 2001;Beckstein and Sansom, 2003;Maibaum and Chandler, 2003). A particular channel conformation may therefore present an energetic barrier to permeation (i.e. be gated closed) in a hydrophobic region of the pore without requiring full steric occlusion. This is referred to as a hydrophobic gate (Aryal et al., 2014b;Beckstein et al., 2001;Jensen et al., 2010;Zhu and Hummer, 2012) or sometimes as a vapour lock (Anishkin and Sukharev, 2004). In such cases, the widening and consequent wetting (i.e. hydration) of the hydrophobic gate region enables the passage of ions through the channel (Fig. 1A).
Computationally, molecular dynamics (MD) simulations of channel structures embedded within a lipid bilayer have aided the functional interpretation of new structures. Such simulations may range in complexity from characterisation of free energy landscapes of ion permeation and/or the conformational transitions associated with gating (Jensen et al., 2010;Zhu and Hummer, 2010) through to simpler simulations of water behaviour within channel pores (Trick et al., 2016).
The influence of pore radius and hydrophobicity on the wetting/de-wetting (i.e. the nanoscale liquid-vapour transition of water inside) of channels have previously been examined by simulation of model nanopores (Beckstein et al., 2001;Beckstein and Sansom, 2003;. For a uniformly hydrophobic constriction, the de-wetted (vapour) state appears stable below a critical pore radius of ~0.5 nm, but this threshold radius decreases if the hydrophobicity of the pore lining is reduced . While channel structures containing hydrophobic regions with radii of ~0.3 nm have been reported to represent de-wetted non-conductive conformations (e.g. (Jia et al., 2018), our ability to predict this based on structure alone is hampered by a lack of detailed understanding of how pore de-wetting depends on the local radius and hydrophobicity in complex biological structures such as ion channels, where a wide range of hydrophobicity profiles are possible via different combinations of pore-lining side chains.
There has been a recent explosion in the structural data available for ion channels, mostly arising from advances in structural biology, especially cryo-electron microscopy (Liao et al., 2014).
Currently there are >800 structures of over 100 unique ion channel proteins deposited in the PDB (Fig. 1B). Furthermore, improved software (e.g. CHAP; www.channotation.org) is now available for analysis of channel pore dimensions, hydrophobicity profiles, and simulated wetting/de-wetting (Klesse et al., 2018). Thus the ion channel structural proteome can now be subjected to a systematic examination of hydrophobic gating. Here we quantify the influence of local pore radius and hydrophobicity on channel de-wetting via MD simulations of water behaviour in nearly 200 ion channel structures. The results provide a structure-based heuristic model enabling rapid prediction of the conductive state of new channel structures almost as soon as they are determined.

A protocol for channel simulation and analysis
Each selected ion channel structure was embedded within a phospholipid bilayer, solvated on either side with 0.15 M NaCl (Fig. 1C), and subjected to three replicate 30 ns atomistic MD simulations to determine the behaviour of water within the transbilayer pore. Positional restraints were applied to backbone atoms to preserve the experimentally determined conformational state of the protein whilst allowing for local side chain flexibility. Using CHAP, side chains that line the pore during simulations were identified and time-averaged profiles calculated for the pore radius, the local hydrophobicity, and the free energy of water as a function of position along the length of the pore ( Fig. 2 and SI Fig. S1). Therefore, each point along the permeation pathway of a channel structure was associated with corresponding values for three variables of interest: pore radius, pore hydrophobicity, and free energy of water at that region of the pore. This enables the dependency of water free energy (and hence pore wetting/de-wetting) on local radius and hydrophobicity to be established and averaged across all the channel structures analysed.

An ion channel dataset
To fully sample the range of radius-hydrophobicity combinations present within the ~810 ion channel structures in the PDB, a reduced dataset was manually curated on the basis of structure quality (primarily resolution) and structural redundancy. Structures with a resolution of 5 Å or worse were excluded, as were those with an incomplete backbone trace in their pore-lining segments. Where several structures of the same ion channel species and of similar pore conformation (judged on root-mean-square-deviations and pore radius profiles along with visual inspection) were available, higher-resolution structures of the wild-type protein were chosen. This yielded a reduced dataset of ~200 structures (see SI Table S1). The MD trajectory dataset (3 x 30 ns for each structure) formed the basis for subsequent analysis of pore de-wetting behaviour.

Analysis of the main energetic barrier in each channel structure
Water free energy profiles derived from equilibrium simulations of wetting/de-wetting within channel structure can be used as the basis for functional annotation (Trick et al., 2016). Any dewetted region presenting an energetic barrier to water is likely to form a barrier to permeation , and thus water permeability may be used as a proxy for ion permeability. Based on our simulation dataset, channel structures containing one or more functionally closed gates were identified. These corresponded to ~70% of the ~200 structures analysed (SI Table S1). In Fig. 3A we show a point for each barrier to water permeability in the (hydrophobicity, radius) plane, coloured by the height of the water free energy barrier (G). Significant barriers (G > 2.5 kJ mol -1 ) are frequently associated with hydrophobic regions. The distribution of water free energy values on the pore hydrophobicity-radius landscape is in qualitative agreement with predictions from simplified models of hydrophobic gating (Beckstein and Sansom, 2003;. Hydrophobic aliphatic side chains (leucine, isoleucine, valine) are frequently found in the pore lining of ion channels at regions with elevated water free energies (Fig. 3B).

Analysis of all residues lining ion channel pores
The relationship between local pore dimensions, hydrophobicity, and free energy of water at any given position along the channel axis becomes clearer when all pore-lining residues in the ~200 structures are surveyed, rather than just those present at barriers in the water free energy profiles.
Examination of the distribution of water density values associated with all pore-lining residues shows a substantive tail in the distribution corresponding to de-wetted regions, i.e. those with reduced water density relative to that of bulk water (SI Fig. S2). If we construct a water free energy landscape in the local pore hydrophobicity-radius plane, a clear pattern emerges (Fig. 4). The landscape is clearly divided into two regions, corresponding to wetted and de-wetted pores. In the hydrophobic region of the landscape an energetic barrier can be encountered in channel regions with radii up to ~0.4 nm (the radius of a water molecule is ~0.15 nm; Fig. 4). Conversely, hydrophilic regions of the landscape include pores which become hydrated at much smaller radii (i.e. < 0.2 nm). We note that the water free energy landscape is very similar even when alternative hydrophobicity scales are employed (SI Fig. S4). By training a support vector machine (SVM) classifier (Cortes and Vapnik, 1995), the hydrophobicity-radius landscape can therefore be divided into two regions corresponding to the likelihood of pore wetting vs. de-wetting.
A heuristic for predicting the conductive state The water free energy landscape derived from the dataset of channel structures and simulations allows us to devise a heuristic technique for predicting the functional state of ion channel structures, in a simulation-free manner based on pore dimensions and hydrophobicity alone (Fig. 5). Having identified pore-lining residues and estimated pore radius and hydrophobicity profiles, the (hydrophobicity, radius) values for the pore-lining residues of a given channel are mapped onto the landscape described above. By identifying the number of residues for which the corresponding (hydrophobicity, radius) points fall below the SVM classification line (i.e. in the de-wetted region of the landscape) one may predict the likelihood of a channel structure corresponding to a de-wetted and hence functionally closed state. The sum of shortest distances for residue points falling below the SVM open vs. closed classification line (Σd) provides a score indicating whether the channel structure is likely to contain a closed gate. As this machine learning-based predictive approach can global_text_v21_figs.docx 16-Dec-18 7 be performed using automated analysis of a single set of coordinates (Klesse et al., 2018) and does not require MD simulation, it can also be performed in a few seconds.
We illustrate this heuristic approach using two recent structures (Fig. 6): one of the TRPV3 channel in a putative sensitised (but non-conductive) conformation (PDB ID: 6MHS) and one of the CRAC channel (Orai) in an open conformation (PDB ID: 6BBF) (Hou et al., 2018). For the open-state Orai channel, none of the (hydrophobicity, radius) points falls below the SVM classification line and so the structure is predicted to be fully wetted and correspond to a functionally open state of the channel. By marked contrast, for the TRPV3 structure, 12 points were below the SVM line, and the sum of shortest distances for residue points falling below this line was Σd = 1.4, leading to clear prediction of a closed hydrophobic gate for this particular channel conformation. Consistent with these predictions, when these two structures were subjected to the complete MD simulation and analysis protocol described above, the resulting water free energy profiles (SI Figure S4) confirmed the heuristic model.

Conclusions
We have used MD simulations of pore water density to determine the wetting/de-wetting properties of nearly 200 ion channel structures. A systematic analysis of the behaviour of pore de-wetting as a function of local radius and hydrophobicity reveals an almost linear dependence of the critical radius for wetting upon the local hydrophobicity. As pore radius and hydrophobicity profiles can now be readily estimated from any known structure, this analysis enables us to propose a simple heuristic model for identifying closed gates in newly determined structures of ion channels (Fig. 7). This model is based on an underlying energy landscape derived from ~600 simulations of water behaviour in ion channels, with a combined duration of ~18 µs. It should be noted that pore wetting of a hydrophobic gate is necessary, but not always sufficient, for ion conduction to occur. However, our method provides a robust exploratory method for the functional annotation of novel channel structures that can be performed prior to more detailed simulations and/or experimental studies of the relationship between ion channel structure and function.
Molecular dynamics simulations were performed with GROMACS (http://www.gromacs.org/) version 5.1 (Abraham et al., 2015), employing the OPLS all-atom protein force field with unitedatom lipids (Jorgensen et al., 1996) and the TIP4P water model (Jorgensen et al., 1983). The integration time-step was 2 fs. A Verlet cut-off scheme was applied, and the Particle Mesh Ewald method (Darden et al., 1993) used to calculate long-range electrostatic interactions. Temperature and pressure were maintained at 37 °C and 1 bar, respectively, using the velocity-rescale thermostat (Bussi et al., 2007) in combination with a semi-isotropic Parrinello and Rahman barostat (Parrinello and Rahman, 1981), with coupling constants of τT = 0.1 ps and τP = 1 ps. Bonds were constrained using the LINCS algorithm (Hess et al., 1997) . A harmonic restraint with a force constant of 1000 kJ mol -1 nm -2 was placed on the protein backbone atoms during simulations. The Channel Annotation Package (CHAP; www.channotation.org; (Klesse et al., 2018)) was used to analyse trajectory frames, with bandwidths of 0.14 nm and 0.45 nm, respectively, applied for estimating water density and hydrophobicity along each channel axis. SVM classification was conducted using the Caret package (Kuhn, 2008) in R version 3.4.4 (www.r-project.org). Figure 1: A Schematic of a hydrophobic gate in an ion channel. A hydrophobic gate region (dashed line) can spontaneously de-wet to form a dry (i.e. vapour) state and functionally close the channel. Widening of the gate leads to wetting (i.e. hydration) of the region and a functionally open state. Water is represented by the pale blue background shading, and ions by red and blue spheres. B Pie chart summarizing the ~200 channel dataset which forms the basis of the current study. The channels are grouped into nine broad families: Cys-loop receptors, ionotropic glutamate receptors, purinergic receptors, transient receptor potential (TRP) channels, K + channels (2TM, 4TM & 6TM), other members of the voltage-gated ion channel (VGIC) superfamily, other cation channels, and other anion channels. C Representation of a typical simulation system used in this study. Shown is the TRPV4 channel in a phospholipid (grey) bilayer. Water molecules and ions are present but omitted for clarity.

Figure 2:
Annotation of an ion channel structure via MD simulations of water within the pore. This is illustrated for the TRPV4 ion channel (PDB ID 6BBJ; (Deng et al., 2018)). A Pore radius profile derived from one of 3 x 30 ns simulations of the protein embedded in a phospholipid bilayer. The mean radius, calculated using the final 20 ns of the trajectory with a sampling interval of 0.5 ns, is shown as a black line (with the grey band representing the standard deviation over time) as a function of position, s, along the pore axis. The pore-lining surface and the pore-lining residues of the channel (hydrophobic in orange; polar in blue) are shown at the top. B Hydrophobicity profile for the pore-lining side chains as a function of position along the pore axis s. The hydrophobicity values are evaluated using a normalised version of the Wimley-White scale (Wimley and White, 1996). C Free energy for water as a function of position along the pore axis s. The free energy profile is obtained from the density of water within the pore as estimated from the MD simulation (see Methods for details).

Figure 3:
A Analysis of the main energetic barrier to water permeation in each ion channel structure. Each point corresponds to a single channel structure (see SI Table S1 for details), indicating the hydrophobicity and pore radius at the highest barrier in the water free energy profile (averaged across the three repeat simulations). The height of the barrier is given by the colour scale. Points with energy values outside the range are shown in the same colour as the lower or upper boundary. For structures containing a narrow pore-loop selectivity filter (e.g. K + channels), the section on their water free energy profile corresponding to the filter region is excluded when locating the maximum energy point to represent the structure. B Frequency of the different pore-lining side chains at an energetic barrier in these ion channel structures, i.e. where the water free energy is greater than 1.5 RT (3.9 kJ mol -1 ) at the mean position of the residue along the pore axis s. Multiple amino acids may occur at each energetic barrier.

Figure 4:
Analysis of all pore-lining residues for all simulated ion channel structures in terms of the local water free energy. From each simulation, mean measurements are noted at the mean position (s) along the pore axis (Fig. 2) of any side chain oriented towards the pore for at least 50% of the simulation. The local water free energy is shown as a function of hydrophobicity and pore radius for all occurrences of pore-lining side chains in the simulated structures. Hydrophobicity values are based on the Wimley-White scale, linearly normalised to [-1, 0.34], in relative units. The hydrophobicity-radius grid is coloured by mean energy. Regions with energy values outside the colour scale range are shown in the same colour as the lower or upper bound. The white dashed contour line indicates the 1 RT (2.6 kJ mol -1 ) position as given by a polynomial support vector machine classifier.

Figure 5:
Schematic of a heuristic prediction approach for the permeation state of an ion channel structure based on analysis of hydrophobicity and pore radius profiles, both of which can be derived from a static structure using the CHAP analysis tool (Klesse et al., 2018).

Figure 6:
Illustration of the heuristic prediction approach. A ROC (receiver operating characteristics) curve analysis based on comparison of simulations of pore wettability (see SI Table S1) and heuristic based predictions indicating an optimal cut-off of Σd = 0.55 for the heuristic prediction of a closed channel, illustrated using two recent structures of (B) the TRPV3 channel in a non-conductive sensitised conformation (PDB ID: 6MHS) and (C) the CRAC channel Orai in an open conformation (PDB ID: 6BBF). Pore radius and hydrophobicity profiles are shown for the transmembrane domains of both structures. For each pore-lining side chain (represented as a coloured square on the pathway profiles), the channel pore radius at the residue is plotted against the corresponding hydrophobicity value. The sum of shortest distances between the dashed (1 RT) contour line and all points falling below it (coloured red), is then 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.

Figure 7:
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 in Fig. 5. Schematic depictions of de-wetted (closed) and hydrated (open) channel are shown for the two main regions of the data. Table S1. Ion channel structures analysed.

Supplementary Information
Detailed list of ~200 channels. Figure S1. Average pore radius, hydrophobicity, and water free energy profiles derived from 3x 30 ns MD simulations of water in the TRPV4 ion channel (PDB ID 6BBJ; Deng et al., 2018).