Cumulative hydropathic topology of a voltage‐gated sodium channel at atomic resolution

Voltage‐gated sodium channels (NavChs) are biological pores that control the flow of sodium ions through the cell membrane. In humans, mutations in genes encoding NavChs can disrupt physiological cellular activity thus leading to a wide spectrum of diseases. Here, we present a topological connection between the functional architecture of a NavAb bacterial channel and accumulation of atomic hydropathicity around its pore. This connection is established via a scaling analysis methodology that elucidates how intrachannel hydropathic density variations translate into hydropathic dipole field configurations along the pore. Our findings suggest the existence of a nonrandom cumulative hydropathic topology that is organized parallel to the membrane surface so that pore's stability, as well as, gating behavior are guaranteed. Given the biophysical significance of the hydropathic effect, our study seeks to provide a computational framework for studying cumulative hydropathic topological properties of NavChs and pore‐forming proteins in general.


Funding information
Seventh Framework Programme, Grant/Award Number: grant n602273; European Union 7th Framework Programme

| INTRODUCTION
Voltage-gated sodium channels (NavChs) are fundamental components of electrically excitable cells such as neurons and muscle cells. 1 They belong to the superfamily of ion channels and their primary function is to facilitate the transport of sodium ions across a cell's membrane in response to changes in the membrane potential. Genetic dysfunction of NavChs has profound implications for physiological neuronal or muscle activity, leading to inherited human disorders. 2 For example, genetic and functional studies have shown that missense mutations in the SCN9A gene encoding the Nav1.7 channel are causally related to neuropathic pain syndromes such as inherited erythromelalgia, [3][4][5][6][7][8][9][10][11][12] paroxysmal extreme pain disorder, [13][14][15][16] and small fiber neuropathy. 17,18 Although molecular mechanisms underlying missense SCN9A mutations remain largely unexplored, structural modeling of the Nav1.7 revealed that disrupting hydropathic lining properties of the Nav1.7's pore can act as a disease-causing molecular trigger. 19,20 The crystal structure of the pre-open I217C NavAb bacterial channel (PDB code: 3RVY) 21  Hydropathic effects not only contribute to the stability of an ion channel's native structure 22,23 but also regulate ion conduction through the pore via "hydropathic gates". [24][25][26][27][28][29] However, the spatial complexity underlying hydropathic interactions 30 Hydropathic moments (usually referred to as "hydrophobic moments" 33 ) provide a simple, yet efficient tool for modeling of molecular mechanisms (eg, see Reference 34). The method of cumulative hydropathic moments 35 has been used for analyzing the spatial profile of cumulative hydropathicity relative to a function-relevant geometry (eg, the geometric center of the residues distribution). [36][37][38][39] Interestingly, it was shown that data extracted from the vanishing, that is, zero-crossing, behavior of the cumulative zero-and secondorder residue hydropathic moments exhibited invariance over a large number of nonredundant globular proteins. 37 This finding highlighted the universality of the spatial transition from a protein's hydrophobic core to its hydrophilic exterior implying the existence of hydropathic scaling laws spanning all different molecular scales.
In principle, hydropathic characteristics along a NavCh's pore are studied at a microscopic scale, that is, by focusing on the hydropathic profile of pore-lining atomic structures (eg, see References 21,40).
In this work, we extend this methodology to incorporate hydropathic characteristics of the spatial transition from the pore's microenvironment toward a macroscopic molecular regime that includes the VSDs. This methodological extension allows for a continuous mapping of hydropathic, as well as, structural properties of a NavCh on two dimensions. Specifically, instead of performing a three-dimensional (3D) mapping of hydropathic characteristics on the NavCh structure at residue-level resolution (as already implemented by available tools, eg, pymol 41 and chimera 42 ), we present a two-dimensional hydropathic mapping procedure of atomic resolution by exploiting poreforming geometrical principles. The main advantage of this approach is that it elucidates how hydropathic dipole accumulation occurs with respect to pore's geometry thus providing information about the pore's gating behavior. 24 In order to achieve that we utilized the tool of cumulative hydropathic moments and developed an atomic sampling algorithm that adapts to geometrical characteristics of the pore so that hydropathic density variations, as well as, their corresponding hydropathic dipole field configurations are visualized and, consequently, analyzed across different molecular scales. Accordingly, the relevant dimensions of the presented methodology are the pore axis coordinate and the molecular scale. We illustrate our methods for the pre-open I217C NavAb channel where we report implications of our observations on molecular stability, as well as, on the formation of hydropathic gates along the pore.
The I217C NavAb model (PDB code: 3RVY) was selected for analysis as it provides a crystallography of the NavAb with the highest resolution (2.7Åresolution). The structure was protonated using the WHAT IF software 43,44 and its principal axes were estimated using the VMD software. 45 A coordinate system b x,b y,b z ð Þwith origin O was introduced and the protonated structure was placed within it, so that the principal pore axis, that is, the axis approximating the direction of the channel's pore, is aligned with the z-axis.

| Geometrical characteristics of the pore
We considered P to represent the principal pore axis and p = (p x , p y , p z ) P to be a pore point. The radius of the smallest sphere that can be squeezed through the pore at p, that is, the pore radius at p, is given by 46 where ‖ Á ‖ is the euclidean norm and vdW i is the Van der Waals radius of the ith atom (see Table S1). Consequently, the distance between p and its nearest neighbor atom is given by In analogy to Equation (1) The unit of measurement for R(p), R p ð Þ , and L(p) is expressed in Å.

| Atomic sampling around the pore
In order to investigate cumulative atomic properties with respect to P, we introduced the p-dependent atomic sampling radius where K α is the total number of sampling spheres centered at p, α denotes the index of the sampling sphere. Note that l α (p) plays the role of the molecular scale, that is, it provides a p-dependent measurement of the size of the channel in Å. The number of sampled atoms within a sphere of radius l α (p) centered at p is given by where θ(Á) is the heaviside function.

| Cumulative hydropathic pore moments
We employed the zero-order hydropathic pore moment. 33 and the first-order (or dipole) hydropathic pore moment 33 where HI χ i = HI i + χ i is the hydropathic value of the ith atom set according to the Kapcha-Rossky atomic hydropathic indices 47 (see Table S2) with additive gaussian noise χ i N μ = 0,σ = 0:001 ð Þ , r ! p,i is the vector from p to c i and the measurement units of h (0) (p, l α (p)) and Þ k are roughly given by kcal/mol and kcal Á Å/mol, respectively. Division of Equations (6) and (7) with Equation (5) provided with an estimation of the hydropathic density 37 and of the hydropathic imbalance 37 pore functions for increasing molecular scale l α (p), respectively. Due to the spherical sampling procedure, all scalar and vector pore functions remain invariant to rotations of the channel around P.

| Discretization of the principal pore axis
We discretized P by introducing the equidistant grid Q = p 1 , 2.6 | Detection of zero-crossing points of a scalar pore function along the principal pore axis We consider f(p, l α (p)) to represent a scalar pore function. If for a given scaling index α, there is a pair {p 0 = p − Δp, p} Q for which the signchange condition f(p 0 , l α (p 0 )) Á f(p, l α (p)) < 0 is satisfied, then the fourdimensional point represents a zero-crossing point of f(p, l α (p)) along p-direction where |Á| returns the absolute value of f and s is obtained by linear interpolation. The set of all detected zero-crossing points of f(p, l α (p)) along p-direction for a given scaling index α is represented as Γ(α). Note that due to the addition of noise to the hydropathic indices, the probability of f(p, l α (p)) being zero was practically eliminated so that f(p

| Hydropathic density variations around preopen I217C NavAb's pore
The first step in this study was to investigate how atomic hydropathicity of the pore' microenvironment changes for increasing molecular scale, l α (p) (see Equation (4)). To do so, we employed the hydropathic density pore function, m (0) (p, l α (p)) (see Equation (8)) and illustrated it as a contour map (see Figure 1A). "Reading" of the contour map was performed from left to right, namely, for a given pore point, p, we investigated the zero-crossing behavior of the hydropathic density pore function for increasing l α (p) based on a color gradient that illustrates negative and positive values of m (0) (p, l α (p)) as blue as red, respectively. Accordingly, blue and red contour map domains represent hydrophobic and hydrophilic atomic structures, respectively. Observations were interpreted with respect to the conserved architectural motif dictating the structural separation of the PDs from the VSDs. 48 For that, a geometrical representation of the structural position of the PDs relatively to the VSDs (and vice-versa) was introduced by approximating the scales ν(p) for which an equilibrium of the corresponding radial distribution functions is achieved (see Suppl.S2).
In that way, the line ν(p) appearing in Figure 1A (and in Figure 2A) indicates that the number of PD atoms within a sampling radius is larger than that of VSD atoms thus roughly approximating the structural boundary between PDs and the VSDs. On the other hand, for l α (p) > ν(p) the number of VSD atoms gradually increases so that the spatial transition from the PDs to the VSDs is realized.
At first glimpse, the contour map of Figure  where sm(ϕ (0) ) represents a smoothed version of ϕ (0) (see Suppl. S4 for construction). Gray shaded areas in A, C mark pore regions where a transition from one pore-lining residue side chain to the next one takes place along P (see Suppl. S5). Dashed black horizontal lines in A and C correspond to extrema of sm(ϕ (0) ). IS stands for intracellular side and ES stands for extracellular side translates into a qualitative change in the profile of the outer pore radius, L(p) (see Equation (3)), as L(p) monotonically decreases and increases for p z ≤ 19.0 and for p z > 19.0, respectively. In that way, a funnel-like outer pore surface with an opening at the IS is formed by 2π-rotation of L(p) around P ( Figure 1A). Deviating from the center of the pore toward the ES, we observe a widening of the pore so that an extracellular opening (or "mouth") is formed (see Figure 1A for p z < −12.0). For l α (p) > 30 Å, hydropathic density variations cannot anymore be deduced via visual investigation of Figure 1A, indicating that the structural combination of the PDs with the VSDs tends to increase hydropathic uniformity of the atomic environment around the pore. This tendency results macroscopically in a down-regulation of hydrophobicity as the channel appears to be weakly hydro- that accounts for the total, that is, whole-channel, hydropathic density.
A detailed examination of hydropathic density variations around the pore was performed by analyzing the zero-crossing behavior of m (0) (p, l α (p)) for increasing molecular scale l α (p). To do so, we detected zero-crossing points of m (0) (p, l α (p)) for increasing scaling index α and constructed the sets Γ (0) (α) (see section 2). A zero-crossing point of m (0) (p, l α (p)) along p-direction identifies a pore point around which the hydropathicity of the atomic environment becomes vanishingly small (see Suppl. S3) and changes polarity, and is illustrated as a blue-to-red contour map transition (see Figure 1A). The arrangement of zerocrossing points on the contour map of Figure 1A revealed the boundaries among visually distinguishable, but also among visually indistinguishable contour domains as indicated by black arrows for l α (p) ≤ ν(p) and l α (p) > ν(p), respectively. Accordingly, the contour map was partitioned into four distinct domains, namely, the hydrophobic contour domain T  Figure 1A). On the other hand, deviating toward the ES or toward the IS we observe an expansion of the PDs, as indicated by the monotonic increase of ν(p) for deviating from the CC, in accordance with an increase in atomic hydrophilicity that is expressed by the formation of Figure 1A). Taken together, these observations reflect the channel's strategy to maintain a hydrophobic core, that is, a hydrophobic CC, by increasing the density of hydrophobic atoms around its mass center, which in turn places hydrophilic atoms elsewhere. The advantage behind this configuration might be an increase in NavAb's stability due to reduction of water-accessible hydrophobic surface achieved by "burying" hydrophobic atoms around the mass center.
In order to obtain a statistical summary of our observations, we introduced the distributions ψ (0) and ϕ (0) , respectively (see Suppl. S4 for construction). This allowed for identifying molecular scales and locations along P where hydropathic density variations tend to increase or decrease as implied by extremization of ψ (0) and ϕ (0) , respectively. As we can see in Figure 1B, ψ (0) can be decomposed into two sub-distributions, namely and ψ The minimum structural information required in order to navigate through the pore's microenvironment is to identify pore-lining residues.
To do that, we partitioned P into distinct pore regions based on a minimal-distance geometrical criterion that takes into account only direct-neighboring residues (see Suppl. S5). This approach revealed that the global minimum of sm(ϕ (0) ) at p z = −12.0 identifies a hydrophilic and narrow pore region that co-localizes with the hydrophilic-to-hydrophilic S178-E177 side-chain transition that is part of the selectivity filter (SF) ( Figure 1A,C and Table S3). The strongly hydrophilic and narrow environment that is locally formed creates favorable conditions for dehydration of incoming ions ( Figure 1A). On the other hand, the global maximum of sm(ϕ (0) ) at p z = 19.0 indicates that closing of the channel at the IS requires a hydropathically diverse atomic environment as indicated by the hydrophilic-to-hydrophobic C217-M221 side-chain transition ( Figure 1A,C and Table S3). The local minimum of sm(ϕ (0) ) at p z = 7.3 identifies the CC and co-localizes with the hydrophobic-to-hydrophobic I210-V 213 side-chain transition ( Figure 1A,C and Table S3). The local maximum of sm(ϕ (0) ) at p z = − 4.4 co-localizes with the hydrophilic-tohydrophobic E177-L176 side-chain transition and identifies a transient pore region located between the narrow, hydrophilic pore region and the CC ( Figure 1A,C and Table S3). Finally, the local maximum of sm(ϕ (0) ) at p z = − 20.3 found at the ES mouth of the pore co-localizes with the hydrophobic-to-hydrophilic M181-S178 side-chain transition ( Figure 1A, C and Table S3). Notably, the hydrophobic-to-hydrophilic M181-S178 side-chain transition becomes evident in Figure 1A as a color change from blue-to-red occurring in the vicinity of the R p ð Þ line with respect to b z ( Figure 1A). Accordingly, the location p z = − 20.3 captures the transition from the wide ES mouth shaped by the hydrophobic M181 side chain toward a narrow and hydrophilic pore region surrounded by the hydrophilic S178 and E177 side chains.

| Cumulative hydropathic topology of pre-open I217C NavAb's pore
In this section, we investigated the dipole field topology underlying hydropathic density variations illustrated in the contour map of Figure 1A. In order to do so, we utilized the hydropathic imbalance pore function, m Þ(see Equation (7)) that can be equivalently understood as the mean (or average) cumulative hydropathic dipole moment (see connection of Equation (7) with Equation (9)) thus providing information about accumulation of hydropathic dipoles around the pore. Given that the gating behavior of a nanopore is largely determined by its geometry and dipole accumulation parallel to its lining, 24 we adopted here a similar approach to the previous section and analyzed how microscopic, pore-lining dipoles scale up and form larger polarized structures that incorporate, both, the PDs and the VSDs. By exploiting radial structural symmetries underlying NavAb's tetrameric conformation it was straight forward to reduce to its z-component, Þ Á b z , as the amplitude of the radial component of dipole hydropathic moment, kh ! xy p, l α p ð Þ ð Þ k(see Equation (7)), is shown to be vanishingly small (see Suppl. S6). Accordingly, our analysis focused on the spatial profile of m Þthat is illustrated on the contour map of Figure 2A. A zero-crossing point of m along p-direction accounts for a change in dipole field's orientation, as well as, for a vanishingly small hydropathic imbalance effect (see Supporting Information S7), that is, for the formation of a hydropathic dipole center on P for a given scaling index α. From a topological viewpoint, this vanishing event corresponds to an equilibrium point of the hydropathic imbalance field acting along P for a given α. Depending on the orientation of the local field, an equilibrium point acts either as molecular attractor (ie, "sink") or molecular repeller (ie, "source") so that a hypothetical molecule (eg, a partially hydrated sodium ion) would be either attracted toward or repelled from it, respectively (see Suppl. S8). The importance of the topological interplay between hydropathic "sinks" and "sources" in ion permeation dynamics and, specifically, ion selectivity was recently highlighted in Reference 49.
Similar to the previous section, we introduced the expression that incorporates topological information extracted from every molecular scale l α (p) where Γ (1) (α) represents the set of all detected dipole centers for a given scaling index α (see section 2). Consequently, the contour map of Figure 2A was partitioned into four domains; the pair T 4 cover the largest contour area so that the map can be roughly split into two parts, namely an ES and an IS part, revealing a topological dichotomy of the pore (Figure 2A). The structural advantage behind this topological configuration is the stabilization of CC via an accumulation of dipoles around it in a direction roughly parallel to P (Figure 2A). Following this line of thought, T  (Figure 2A,B), and the topological dichotomy of the pore is disrupted. Specifically, the structural combination of the PDs with the VSDs causes a displacement of dipole centers toward the IS, as indicated in Figure 2A by the "(*)" arrow, so that the size of T  Table S3). This topological configuration facilitates a passage for an incoming hydrated ion through the ES mouth toward the SF under the influence of a strong, bi-directional field effect originating from the combination of T 3 . The local maximum of sm(ϕ (1) ) at p z = − 11.7 co-localizes with the narrowing of the pore thus capturing the local minimization of R(p), as well as, with the global minimum of sm(ϕ (0) ) ( Figures 1A and 2A,C). Given that dipole accumulation around p z = − 11.7 occurs roughly parallel to P, the relative orientation of T  Table S2) performed in pymol 41 computational environment. B, Conservation analysis of a PD structural unit. Analysis was performed online on the ConSurf Server (http://consurf.tau.ac.il/) according to the algorithmic implementations described in References 55,56. C, A summary read-out of cumulative hydropathic topological analysis is provided in terms of the strongly and weakly smoothed topological quasiprobability score (TQS) (see Suppl. S8, Equation (S12)) plotted for p Q. Geometrical pore characteristics are represented in terms of the pore radius trace, R(p), plotted for p Q. Pore points belonging to the subsets Q 2Ái + 1, i = 0, 1, …, 7 minimize their distance from the pore-lining residue side chains M181, S178, E177, L176, I210, V 213, C217 and M221, respectively. Vertical green zones in C highlight the locations of persistent attractors and repellers along P corresponding to negative-value and positive-value maxima of TQS, respectively (see also Suppl. S8), and indicating the location and relative size of energy barriers imposed to permeating molecules (see Suppl. S8). ES stands for extracellular side, EF for extracellular funnel, SF for selectivity filter, CC for central cavity and AG for activation gate. Note that the SF residue side chains complex T175-L176-E177-S178-W179 averages the highest conservation score with 2.4 (out of 4.0) with porelining residue conservation scores confined within [0, 2.0] [Color figure can be viewed at wileyonlinelibrary.com] the spatial organization of the SF T175-L176-E177-S178-W179 residue complex 21 and, especially, of E177 side chains playing a key role in ion selectivity in sodium and calcium channels. 54 The local maximum of sm(ϕ (1) ) at p z = 22.7 identifies with high precision the pore point where pore occlusion takes place revealing a microscopic dipole accumulation mechanism occurring roughly parallel to P, as we can visually deduce from Figure  What are the lessons to be learned from cumulative hydropathic topological analysis? A summarizing read-out of the cumulative hydropathic topological analysis is presented in Figure 3 alongside with targeted conservation analysis, 55,56 as well as, 3D hydropathicity mapping of the pore walls. Specifically, a condensed output of our method was obtained in terms of a topological quasiprobability score (TQS) (see Suppl. S8, Equation (S12)) that quantifies the probability along the pore for a molecular attractor (or repeller) to occur ( Figure 3C).
TQS essentially provides with a qualitative description of the relative size and location of energetic barriers to molecular permeation, that is, "gates," 29 along the pore (see Suppl. S8). According to this interpretation scheme, a central gating mechanism keeps the hydrophobic CC "closed" by imposing large energy barriers to sodium ions entering from the ES or the IS side while facilitating an "opening" to small hydrophobes and waters entering via the side fenestrations. Hence, hydrated sodium ions escaping from the SF' local binding sites have to overcome a large energy barrier in order to arrive at the center of the pore ( Figure 3C). This barrier stems from the molecular influx into CC through the fenestrations which is in opposite direction to the ion escape trajectory thus preventing (or, at least, not favoring) ion entrance into the CC. Escape trajectory is thus likely to involve multiple, partial re-and dehydration cycles, that is, an interplay of attractors-vs-repellers sites, so that ions are gradually transported from the narrow, strongly hydrophilic and highly conserved SF environment toward the wide, hydrophobic center of the pore (see The advantage behind a long escape trajectory is that NavAb's specificity is optimized as entrance into CC is permitted only to ions that can re-and de-hydrate according to the local gradients of the energy landscape. Ions arriving in the center of the pore are fully hydrated thus they can rapidly diffuse toward the hydropathically diverse AG where a local attractor site is formed at p z ≈ 19.0 (see Figure 3A,C where TQS attains its local negative maximum, and previous Section for properties of the atomic environment around p z ≈ 19.0).
Finally, a weak, repeller mechanism at p z ≈ − 20.3 attributed to the hydrophobic-to-hydrophilic interplay among the M121 and S178 side chains is likely to regulate ion access to the SF from the ES (see Figure 3A,C where TQS attains its local positive maximum, and previous Section for properties of the atomic environment around p z ≈ − 20.3).

| DISCUSSION
In this study, we followed the trace of earlier works on cumulative hydropathic analysis of protein systems (see 50) and presented a computational scheme that allows for visual detection and analysis of hydropathic pattern formation around a NavCh's pore. Specifically, considering that structural symmetries play a fundamental role in biological pore (for a review, see Reference 57), we developed a scaling analysis methodology that utilizes the tools of hydropathic density and hydropathic imbalance pore functions (see Equations (8) and (9), respectively) and returns two-dimensional spatial profiles (ie, contour maps) of these experimental quantities where the relevant dimensions are unfolding perpendicular and parallel to the membrane surface. In the absence of structural pore symmetries, for example, in the case of heteromeric eukaryotic NavChs such as the Nav1.7, our algorithm is still applicable but radial contributions to the hydropathic dipole field need to be taken into account. This can be done by introducing an appropriate color scheme for illustrating changes in orientation and magnitude of the radial hydropathic dipole field component (see Equation (7)) so that a third contour map is con- We demonstrated that the spatial profiles of hydropathic density and hydropathic imbalance pore functions exhibit pseudo-symmetrical characteristics with respect to an axis parallel to the membrane surface so that the pore is roughly dichotomized with respect to the CC. This becomes evident if we focus on the relative positioning of the pairs T . The biophysical principle underlying this mode of organization is the optimization strategy dictating that burying of hydrophobic atoms around a protein's core contributes to protein stability 58 (the pre-open I217C NavAb's hydrophobic core was identified as the CC (see Figure 1A)).
Persistent extrema of the distributions ϕ (0) and ϕ (1) were utilized in order to investigate hydropathic gating properties of the pore. In particular, we demonstrated that the maxima of sm(ϕ (1) ) identify pore locations where hydropathic topological and geometrical characteristics of the atomic environment favor the formation of a central hydrophobic gate (at p z ≈ 2.8) regulating ion transport between the hydrophilic SF attractor (at p z ≈ − 11.7) and the hydropathically diverse activation gate (AG) 21 (at p z ≈ 22.7). Conservation analysis of the PDs highlighted the highly conserved nature of the SF as already pointed out in Reference 21 (see Figure 3). Accordingly, the SF attractor topological configuration reflects an evolutionary-driven functional optimization of the PDs that might be relevant for other NavCh species as well. Given that these locations account for a smoothed, multiscale hydropathic effect, thermal fluctuations are expected to induce only small dislocations.
In summary, the strength of the presented methodologies is that they allow not only for a detailed mapping of multiscale hydropathic characteristics lining a NavCh's pore but also for elucidating how their underlying topology correlates with the pore's geometry at atomic resolution. This is a crucial step toward understanding where along a NavCh's pore hydropathic gates are likely to form and how they contribute to its functioning. A weaknesses, however, is that multiscale hydropathic mapping at atomic resolution is computationally expensive in comparison to traditional structural biology methods that operate at residue level and that interpretation of observables can be far from trivial. To circumvent these complexities statistical summary measures, such as the distributions ψ and ϕ, can be employed alongside with targeted conservation analysis. Xenakis wrote the manuscript in consultation with all the co-authors; PROPANE study group provided with technical and scientific support.

ACKNOWLEDGMENTS
The study was partly funded by the European Union 7th Framework Programme (grant n602273).