Cumulative hydrophobic topology of the NavAb's pore ...

Cumulative hydrophobic topology of the NavAb’s pore at atomic resolution. Xenakis M.N.1,2,* , Kapetis D.3 , Yang Y.4,5 , Heijman J.6 , Waxman S.G.7,8 , Lauria G.3,9 , Faber C.G.10 , Smeets H.J.M.2 , Westra R.L.1 , Lindsey P.2 and PROPANE Study Group** 1

Department of Data Science and Knowledge Engineering, Maastricht University, PO Box 616, 6200 MD Maastricht, the Netherlands 2 Department of Clinical Genetics, Maastricht University, PO Box 616, 6200 MD Maastricht, The Netherlands 3 Unit of Neuroalgology, IRCCS Foundation ”Carlo Besta” Neurological Institute, via Celoria 11, 20133 Milan, Italy 4 Department of Medicinal Chemistry and Molecular Pharmacology, Purdue University College of Pharmacy, West Lafayette, IN, 47907, USA 5 Purdue Institute for Integrative Neuroscience, West Lafayette, IN 47907, USA 6 Department of Cardiology, CARIM School for Cardiovascular Diseases, Maastricht University, PO Box 616, 6200 MD Maastricht, The Netherlands 7 Department of Neurology and Center for Neuroscience and Regeneration Research, Yale University School of Medicine, New Haven, CT 06510, USA. 8 Rehabilitation Research Center, Veterans Affairs Connecticut Healthcare System, West Haven, CT 06516, USA. 9 Department of Biomedical and Clinical Sciences ”Luigi Sacco”, University of Milan, via G.B. Grassi 74, 20157 Milan, Italy 10 Department of Neurology, Maastricht University Medical Center, PO Box 5800, 6202 AZ Maastricht, The Netherlands

June 3, 2018

Abstract 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 disturb heart rhythm and affect the nervous system. NavChs are pharmacological targets of numerous drugs, thus elucidating their structure and function is of major importance for developing safer, more * Correspondence and requests for materials should be addressed to M.N.X. E-mail: [email protected] ** PROPANE Study Group (in alphabetical order): Al Momani R., Boneschi F.M., Cazzato D., Cestle S., Chever O., Clarelli F., Dib-Hajj S.D., Eijkenboom I., Gerrits M.M., de Greef B., Hoeijmakers J.G.J., Lombardi R., Lopez I., Mailk R., Mantegazza M., Marchi M., Merkies I.S.J., Quattrini A., Santoro S., Sopacua M., Szklarczyk R., Taiana M, Tavakoli M., Vanoevelen J.M., Zauli A., Ziegler D.

1

effective modulators. Here, we present a topological connection between the functional architecture of the bacterial voltage-gated sodium channel NavAb crystallized at the pre-open state and the spatial organization of hydrophobicity around its pore at the atomic scale. Our analysis reveals the existence of a non-random cumulative hydrophobic topology that is organized parallel to the membrane. This constitutes a novel type of emergent property of the NavAb channel and encodes information about the functional modularity of the pore domains, as well as, about hydrophobic stability effects mediated by the voltage sensing domains. Given the fundamental biophysical importance of hydrophobic interactions, our algorithm seeks to provide a basis for novel in silico tools for structural analysis of NavChs and pore-forming proteins in general.

1

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 regulate the transport of sodium ions across a cell’s membrane in response to changes in the membrane potential. Dysfunction of NavChs due to mutations in the underlying genes has profound implications for physiological neuronal or muscle activity, leading to inherited human disorders such as chronic pain, cardiac arrhythmias and epilepsy [2]. For example, genetic and functional studies have shown that gain-of-function mutations in the SCN9A gene encoding the Nav1.7 channel are causally related to neuropathic pain syndromes such as inherited erythromelalgia [3–12], paroxysmal extreme pain disorder [13–16] and small fibre neuropathy [17, 18]. The crystal structure of the NavAb at the pre-open state [19] provides a prototype for understanding the functional architecture of the NavChs family. The NavAb structural model consists of four homologous subunits, each of which comprises a pore domain (PD) interlinked with a voltage-sensing domain (VSD). Interestingly, an interplay of hydrophobic and hydrophilic residues results in a pore structure where functionally distinct regions coexist [19]. This finding highlights the fundamental role that hydrophobic effects play in shaping the NavAb’s functional architecture but also raises the question how this ensemble of heterogeneous structural components stabilizes within the membrane. Hydrophobic effects not only contribute to the stability of an ion channel’s native structure [20,21] but also regulate ion conduction through its pore [22–25]. In humans, gain-of-function mutations in the PDs can perturb atomic hydrophobic channel properties and, consequently, alter channel gating so that disease symptoms are triggered. For instance, disruption of a hydrophobic ring formed by four pore-lining aromatic residues at the putative activation gate (AG) of the Nav1.7 has been proposed as a plausible disease-causing mechanism associated with inherited forms of neuropathic pain [26, 27]. Furthermore, hydrophobic effects are crucial for drug activity as they are known to regulate drug-channel

2

interactions in the pore [28]. Although there is a tremendous interest in obtaining a better understanding of the role of hydrophobic effects in NavChs under pathophysiological conditions, the spatial complexity underlying hydrophobic interactions poses a significant obstacle in the development of computational tools for molecular hydrophobicity analysis [29–31]. Hydrophobic moments [32] provide a simple but efficient tool for modeling of complex molecular mechanisms (e.g., see [33]). In principle, hydrophobic characteristics along a NavCh’s pore are studied at a microscopic scale, i.e., by focusing on atomic structures lining the pore (e.g., see [19, 34]). Although this approach is useful for studying short-range hydrophobic interactions, it disregards longrange hydrophobic contributions to the pore microenvironment originating from spatially-distant components such as the voltage-sensing domains. In this work, we attempt to bridge this spatial gap by presenting a methodology for analysis of hydrophobic pattern formation around a NavCh’s pore. To do so, we utilized the method of cumulative hydrophobic moments [35–38] and developed a nested atomic sampling algorithm that adapts to geometrical characteristics of a NavCh’s pore so that a spatially-adjusted computation of cumulative atomic properties is performed. As we demonstrate in this report, extending this algorithmic procedure by a topological analysis step, allows for extracting stabilityand function-relevant channel information. This is illustrated for the pre-open NavAb channel where cumulative hydrophobic characteristics were quantified in terms of the zero- and first-order cumulative hydrophobic moment functions.

2 2.1

Results Accumulation of hydrophobicity around NavAb’s pore

The first step in this study was to investigate how, i.e., in what kind of spatial pattern, atomic hydrophobicity accumulates around NavAb’s pore. To do so, we employed the zero-order mean cumulative hydrophobic pore moment, m(0) (p, lα (p)) (see equation m8), and analyzed its spatial profile that is illustrated as a contour map in Fig. 1(a). Cumulative hydrophobic patterns observed in Fig. 1(a) were interpreted with respect to the structural organization of the PDs and VSDs around P . For that, a geometrical representation of the relative positioning of the PDs with respect to the VSDs (and vice-versa) was introduced by approximating the distance ν(p) from p at which the structural transition from the PDs to the VSDs takes place (SI S2). Accordingly, the line ν(p) appearing in Fig. 1(a) indicates that the number of PD atoms within a sampling sphere of radius lα (p) < ν(p) is always larger than that of VSD atoms. On the other hand, for lα (p) > ν(p) the number of VSD atoms gradually increases so that a mixture of PD and VSD atoms is obtained. At first glimpse, the contour map of Fig. 1(a) reveals that the spatial organization of PD atoms around P results in the formation of a centrally-located hydrophobic cavity. This becomes evident if we focus on the contour area enclosed within −5.3 ≤ pz ≤ 6.6 and lα (p) < 13.5 ˚ A (Fig. 1(a)). Deviating from this

3

area toward the IS, we observe the narrowing and, eventually, the geometrical closing of the pore by a hydrophilic pore-lining structure as we can see in Fig. 1(a) for pz > 18.0. Noticeably, the geometrical closing of the channel appears to correlate with the profile of the outer pore radius as L(p) monotonically decreases for pz ≤ 18.0, while for pz > 18.0 the situation is reversed as L(p) increases toward the IS (Fig. 1(a)). Accordingly, a smooth, funnel-line outer pore surface is obtained by 2π-rotation of L(p) around P , which is characterized by an opening coinciding with the closing of the pore (Fig. 1(a)). At the ES, we observe a widening of the pore that is lined by mainly hydrophobic atoms as shown in Fig. 1(a) for pz < −20.1. Obviously, for lα (p) = L(p) the contour map of Fig. 1(a) has become uniform as individual traces of m(0) (p, lα (p)) converge to (0) mtotal =m(0) (p, L(p)) =Nc−1

Nc X

HIξ,i = 0.00487 ∼ kcal/mol

(r1)

i

that represents the total, i.e., whole-channel, mean hydrophobicity. Noticeably, uniformization of the contour map coincides with the structural transition from the PDs to the VSDs as indicated by dissolution of cumulative hydrophobic patterns for lα (p) ≥ ν(p) (Fig. 1(a)). A detailed examination of contour pattern formation was performed by analyzing the zero-crossing behavior of m(0) (p, lα (p)) across different molecular scales. To do so, we detected zero-crossing points of m(0) (p, lα (p)) for every α ∈ A and constructed the sets Γ(0) (α) (see Methods). A zero-crossing point of m(0) (p, lα (p)) along p-direction is represented as (s(0) , lα (s(0) )) and identifies a location on P around which hydrophobicity originating from an atomic cluster of size lα (s(0) ) becomes vanishingly small (SI S3). Thus, it accounts for a transition from a hydrophobic to a hydrophilic, i.e., hydrophobic-to-hydrophilic, (or from a hydrophilic to a hydrophobic, i.e., hydrophilic-to-hydrophobic) cluster of atoms along P , which is visualized as a blue-to-red (or red-to-blue) color change in Fig. 1(a). Consequently, the mutual arrangement, i.e., the topology, of zero-crossing points of m(0) (p, lα (p)) encodes information about formation and dissolution of contour patterns. This motivated us to introduce the union of all Γ(0) (α) sets [ Ω(0) = {Γ(0) (α)} (r2) α

that serves as a topological representation of the spatial profile of m(0) (p, lα (p)).

4

Figure 1: Spatial profile of the zero-order mean cumulative hydrophobic pore moment. (a), Contour map of m(0) (p, lα (p)) for p ∈ Q and α ∈ A. Blue and red color contour domains represent hydrophobic and hydrophilic cu¯ mulative atomic structures around P , respectively. Black lines R(p), R(p) and L(p) depict geometrical pore characteristics. Magenta dashed line ν(p) serves as a geometrical representation of the structural transition from the PDs to the VSDs (SI S2). Black arrows indicate the clustering behavior of zero-crossing points of m(0) (p, lα (p)) found in Ω(0) resulting in the formation of the contour (0) (0) (0) (0) domains T1 , T2 , T3 and T4 . (b), Histogram of Ω(0) along lα (p)-direction, ψ (0) (SI S5). (c), Histogram of Ω(0) along p-direction, φ(0) (SI S5). sm(φ(0) ) represents a smoothed version of φ(0) (SI S5). M 181, S178, E177, L176, I210, V 213, C217, M 221 represent pore-lining residue side chains (SI S4). Grey shaded areas in (a),(c) mark pore regions, where a structural transition from one pore-lining residue side chain to the next one takes place along P (SI S4). Dashed black vertical and horizontal lines in (a),(c) correspond to local extrema of sm(ψ (0) ) and sm(φ(0) ), respectively. IS stands for intracellular side and ES stands for extracellular side. As we show in Fig. 1(a), plotting of the zero-crossing points found in Ω(0) highlights 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 of Fig. 1(a) can be partitioned into four topologically distinct domains, namely,

5

(0)

(0)

the hydrophobic contour domain T3 , the hydrophilic contour domains T1 (0) (0) and T2 , and the weakly hydrophilic contour domain T4 . This topological treatment allowed for deducing spatial correlations that would otherwise remain elusive. In particular, it became obvious that the centrally-located hydrophobic cavity is a smaller part of a larger domain that covers the largest area of the (0) (0) (0) (0) contour map, namely of T3 (Fig. 1(a)). Around T3 the hydrophilic T1 , T2 (0) and T4 are formed in accordance with the overall structural organization of the PDs and the VSDs. In particular, the spatial acceleration of the structural transition from the PDs to the VSDs, i.e., minimization of ν(p), takes place for −5.3 ≤ pz ≤ 6.6 so that it coincides with the formation of the centrally-located hydrophobic cavity (Fig. 1(a)). On the other hand, deviating toward the ES or toward the IS we observe the spatial deceleration of the structural transition from the PDs to the VSDs, i.e., ν(p) monotonically increases, in accordance with an increase in cumulative hydrophilicity that is topologically expressed by (0) (0) the formation of T1 and T2 (Fig. 1(a)). Taken together, these observations reflect the channel’s tendency to maintain a hydrophobic core by increasing the density of hydrophobic atoms around its mass center, which in turn places hydrophilic atoms elsewhere. This resembles a stability mechanism where a local hydrophobic inhomogeneity is wiped out so that the relative stability of spatially-distant structural components is ensured. Accordingly, the succession (0) (0) of T3 by T4 in Fig. 1(a) can be understood as a counteract to the preferential packing of hydrophobic atoms around the channel’s center, which increases whole-channel stability via an up-regulation of total hydrophilicity as indicated (0) by mtotal . In order to summarize topological information along lα (p)- and p-direction, we introduced the topological distributions ψ (0) and φ(0) , respectively (SI S5). This allowed for identifying molecular scales and locations on P where hydrophobic variability, i.e., the number of zero-crossing points, tends to increase or decrease in terms of ψ (0) and φ(0) , respectively. As we can see in Fig. 1(b), ψ (0) can be decomposed into two sub-distributions encoding qualitatively different (0) (0) topological characteristics, namely, sub-distributions ψA and ψB (Fig. 1(b)). (0) (0) A comparison of ψA with ψB reveals that the hydrophobic variability within the PDs is significantly increased as the majority of zero-crossing points (approximately 65%) are found in the contour map for lα (p) < ν(p). The coexistence (0) (0) of ψA with ψB reflects the distinctive structural and hydrophobic roles of the PDs and of the VSDs and, consequently, identifies the interface between them, i.e, the set of scales characterized by lα (p) ≈ ν(p), as a structural barrier where hydrophobic variability disappears (Fig. 1(b),(a)). Given that pore-lining residues play a key role in shaping the physicochemical properties of the NavAb’s pore microenvironment [19], we introduced a partitioning of the pore based on a minimal-distance geometrical criterion that mapped pore-lining residues on it (SI S4). This enabled for interpreting our findings with respect to the structural sequence of pore-lining residues along while focusing on locations on P where multiscale hydrophobic variability is

6

maximized (or minimized) in terms of the smoothed profile of φ(0) , sm(φ(0) ) (SI S5). We found that the global minimum of sm(φ(0) ) at pz = −12.5 identifies a hydrophilic and narrow pore region that co-localizes with the hydrophilic-tohydrophilic S178-E177 transition (Fig. 1(a),(c) and SI Table S3). This finding (0) indicates that T1 initiates as a microscopic, predominantly hydrophilic atomic cluster formed at the structural interface between the hydrophilic side chains S178 and E177 and gradually expands until it reaches a critical cluster size, (0) (0) followed by T3 and T4 . On the other hand, the global maximum of sm(φ(0) ) (0) at pz = 18.0 accounts mainly for the formation of T2 and co-localizes with the hydrophilic-to-hydrophobic C217-M 221 transition (Fig. 1(a),(c) and SI Table ¯ S3). Interestingly, as we can see in Fig. 1(a) in the vicinity of R(p) line, hydrophilic and hydrophobic atomic components of C217 and M 221 side chains, respectively, are preferentially placed in the vicinity of P so that the geometrical closing of the channel is achieved. The local minimum of sm(φ(0) ) at pz = 6.6 identifies the centrally-located hydrophobic cavity and co-localizes with the hydrophobic-to-hydrophobic I210-V 213 transition (Fig. 1(a),(c) and SI Ta(0) ble S3). Accordingly, the microscopic origins of T3 can be traced down to the structural alignment of the I210 and V 213 side chains along P . The local maximum of sm(φ(0) ) at pz = −5.3 accounts mainly for the topological bound(0) (0) ary between T1 and T3 and co-localized with the hydrophilic-to-hydrophobic S178-E177 transition. Thus, it identifies a transient pore region located between the narrow, hydrophilic pore region and the centrally-located hydrophobic cavity (Fig. 1(a),(c) and SI Table S3). The local maximum of sm(φ(0) ) at pz = −20.1 occurs at the ES mouth of the pore and co-localizes with the hydrophobic-tohydrophilic M 181-S178 transition (Fig. 1(a),(c) and SI Table S3). It results (0) (0) from the topological boundary between T1 and T3 , as well as, from irreg(0) (0) ular clustering of zero-crossing points between T3 and T4 . Noticeably, the hydrophobic-to-hydrophilic M 181-S178 transition becomes evident in Fig. 1(a) ¯ as a color change from blue-to-red occurring in the vicinity of the R(p) line with respect to ˆ z (Fig. 1(a)). These events reflect the transition from a wide and hydrophobic ES channel’s mouth shaped by the hydrophobic M 181 side chains to a narrow and hydrophilic pore region surrounded by the hydrophilic S178 and E177 side chains.

2.2

Accumulation of hydrophobic dipoles on NavAb’s pore

Next, we investigated hydrophobic dipole formation and accumulation around P ~ (1) (p, lα (p)) in terms of the first-order mean cumulative hydrophobic pore moment, m (1) ~ (p, lα (p))|| (see equation m9). From a structural biology point of view, ||m can be understood as a measure of the hydrophobic imbalance [38] of the surrounding atomic structures with respect to p. Due to radial structural symmetries underlying the channel’s tetrameric conformation, we were able to reduce (1) (1) ~ (1) (p, lα (p)) to its z-component m ~ z (p, lα (p)) = mz (p, lα (p)) · ˆ m z (SI S6). (1) Accordingly, our analysis focused on the spatial profile of mz (p, lα (p)), which

7

accounts for multiscale hydrophobic imbalance effects acting parallel to P (Fig. 2(a)). In order to extract biophysically-relevant information from the contour map of Fig. 2(a) we followed the same steps as in the previous section and in(1) vestigated the zero-crossing behavior of mz (p, lα (p)) by constructing the sets (1) Γ(1) (α) (see Methods). A zero-crossing point of mz (p, lα (p)) along p-direction (1) (1) is represented as (s , lα (s )) and accounts for a canceling out of the hydrophobic imbalance (SI S7), i.e., for the formation of a hydrophobic dipole centered [39] at s(1) . Consequently, s(1) identifies a location on P where an equilibrium or ”balanced” hydrophobic configuration of the surrounding atomic cluster of size lα (s(1) ) is achieved. Dipole centers are illustrated in the contour map of Fig. 2(a) as a blue-to-red (or red-to-blue) color changes and, from a topological point of view, represent a source or a sink of the vector field ~ (1) (p, lα (p)) along p-direction, i.e., a change in its orientation from pointm ing toward the ES to pointing toward the IS or from pointing toward the IS to pointing toward the ES, respectively. Similarly to the previous section, we investigated spatial patterns of dipole formation in terms of topological properties of the cumulative hydrophobic dipole field, which were summarized with the expression [ Ω(1) = {Γ(1) (α)} (r3) α (1)

As shown in Fig. 2(a) the spatial profile of mz (p, lα (p)) is underlined by a complex, but nevertheless non-random, cumulative hydrophobic topology that is organized parallel to the membrane. In particular, the contour map of Fig. 2(a) (1) (1) can be roughly partitioned into four contour domains, namely, the T1 , T4 and (1) (1) (1) ~ z (p, lα (p)) T2 , T3 contour domains that correspond to configurations of m (1) with orientation toward the ES, i.e., mz (p, lα (p)) < 0, and with orientation (1) (1) (1) toward the IS, i.e., mz (p, lα (p)) > 0, respectively. Strikingly, T3 and T4 cover the largest area of the contour map and roughly split it into two an ES and an IS part (Fig. 2(a)). This topological configuration reveals that the increased density in hydrophobic atoms around the centrally-located cavity is stabilized by an accumulation of dipoles taking place roughly perpendicular to P , i.e., parallel to the membrane, and expanding up to lα (p) ≈ ν(p) (Fig. 2(a)). Following this (1) (1) line of thought, the formation of T1 and T2 accounts for stabilization of hydrophilic atomic structures at the ES and IS associated with the formation (0) (0) of T1 and T2 , respectively (Fig. 1(a), 2(a)). Interestingly, in contrast to (0) (1) ψ , ψ forms a single distribution indicating that the vast majority of dipole centers (approximately 90%) are characterized by lα (s(1) ) < ν(p) highlighting the prominent role that PDs play in shaping the dipole force field acting along P (Fig. 2(b)). For lα (p) ≈ ν(p) a qualitative change in the spatial profile of (1) mz (p, lα (p)) takes place where the structural transition from the PDs to the VSDs is accompanied by a decrease in the number of dipole centers as indicated by flattening of ψ (1) (Fig. 2(a),(b)). Specifically, the structural combination of the PDs with the VSDs causes a displacement of dipole centers toward the 8

IS, as indicated in Fig. 2(a) by the ”(d)” arrow, reflecting the up-regulation of total hydrophilicity (see equation r1). For lα (p) = L(p) a unidirectional field configuration has emerged with (0)

m(1) z (p, L(p)) = −mtotal · pz + β

(r4)

that is in accordance with the molecular moments field theory presented in [39] (SI S8).

Figure 2: Spatial profile of the first-order mean cumulative hydropho(1) bic pore moment. (a), Contour map of mz (p, lα (p)) for p ∈ Q and α ∈ A. (1) ~ z (p, lα (p)) with Blue and red contour domains represent configurations of m orientation toward the ES and toward the IS, respectively. Black lines R(p), ¯ R(p) and L(p) depict geometrical pore characteristics. Magenta dashed line ν(p) serves as a geometrical representation of the structural transition from the PDs to the VSDs (SI S2). Black arrows indicate the clustering behavior of (1) dipole centers found in Ω(1) resulting in the formation of contour domains T1 , (1) (1) (1) T2 , T3 and T4 . Black arrow labeled as ”(d)” indicates the clustering behavior of dipole centers for lα (p) > ν(p). (b), Distribution of hydrophobic centers found in Ω(1) along lα (p)-direction, ψ (1) (SI S5). (c), Distribution of hydrophobic centers found in Ω(1) along p-direction, φ(1) (SI S5). sm(φ(1) ) represents a smoothed version of φ(1) (SI S5). M 181, S178, E177, L176, I210, V 213, C217, M 221 represent pore-lining residue side chains (SI S4). Grey shaded areas in 9

(a),(c) mark pore regions, where a structural transition from one pore-lining residue side chain to the next one takes place (SI S4). Vertical and horizontal dashed black lines in (a),(b),(c) correspond to local extrema of the smoothed distribution sm(ψ (1) ) and sm(φ(1) ), respectively. IS stands for intracellular side and ES stands for extracellular side. Analogous to the previous section, we employed the smoothed profile of the topological distribution φ(1) , sm(φ(1) ) (SI S5), in order to identify key locations on P that merge multiscale structural and hydrophobic information. Interestingly, we found that the global maximum of sm(φ(1) ) at pz = 2.7 not only co-localizes with the hydrophobic-to-hydrophobic L176-I210 transition but also with the local maximization of R(p), both identifying the centrally-located hydrophobic cavity (Fig. 2(a),(c) and SI Table S3). This finding demonstrates how microscopic channel characteristics scale up and form larger structures. Specif(0) ically, T3 appears to be stabilized by an accumulation of dipoles that initiates around the pore region lined by the hydrophobic side chains L176 and I210 and (1) (1) gradually expands so that the pair {T3 , T4 } is formed (Fig. 1(a), 2(a)). On the other hand, the global minimum of sm(φ(1) ) at pz = −17.9 accounts for the (1) formation of T1 and co-localizes with the hydrophobic-to-hydrophilic M 181S178 transition that signals the narrowing of the pore (Fig. 2(a),(c) and SI Table S3). The local maximum of sm(φ(1) ) at pz = −11.7 accounts for the topo(1) (1) logical boundary between T1 and T3 and it almost overlaps with the global minimum of sm(φ(0) ) (Fig. 1(a), 2(a),(c) and SI Table S3). This implies that (0) T1 is stabilized by a dipole accumulation initiating at the interface between the hydrophilic side chains S178 and E177 and expanding up to lα (p) ≈ 13.5 ˚ A. The (1) (1) local maximum of sm(φ ) at pz = 22.7 accounts for the formation of T2 but (1) (1) also for the topological boundary between T3 and T4 and co-localizes with the hydrophobic M 221 side chain (Fig. 2(a),(c) and SI Table S3). Strikingly, it coincides with the geometrical closing of the channel where an accumulation of dipoles perpendicular to P is observed as we can see in Fig. 2(a) in the vicinity ¯ of the R(p) line. This microscopic event appears to stabilize the alignment of hydrophilic M 221 atoms in the vicinity of P so that closing of the channel is achieved (Fig. 2(a)). The local minima at pz = −7.1 and pz = 10.8 identify locations on P where the orientation of the hydrophobic imbalance vectors remains (1) (1) fixed due to formation of T3 and T4 and co-localize with the hydrophilic-tohydrophobic E177-L176 and with the hydrophobic-to-hydrophilic V 213-C217 transitions, respectively. Accordingly, these locations merge information about the transition from the hydrophobic and wide pore center toward pore regions of qualitatively different structural and hydrophobic characteristics (Fig. 2(a),(c)).

3

Discussion

In this study we followed the trace of the early works of Silverman (see [36–39]) and presented a systematic methodology that allows for investigating cumula10

tive hydrophobic pattern formation around a channel’s pore. Specifically, we utilized the concepts of hydrophobicity density [38] and hydrophobic imbalance [36] in order to obtain estimations of cumulative hydrophobic effects that are independent of atomic density variations. Thus, reported contour patterns account solely for cumulative hydrophobic effects emerging from the spatial organization of atomic structures around P . For the sake of simplicity and due to limitations in computing time, we treated NavAb as a rigid entity by neglecting thermal atomic fluctuations. Nevertheless, introduction of thermal fluctuations is expected to perturb microscopic hydrophobic characteristics, but not affect qualitative, i.e., macroscopic, contour pattern features such as the emergence of topological contour domains. (1) We demonstrated that m(0) (p, lα (p)) and mz (p, lα (p)) are underlined by topologies that exhibit pseudo-symmetrical characteristics with respect to an axis parallel to the membrane roughly dichotomizing the pore (note that the word ”pseudo” is used here as a synonym to imperfect). This becomes evi(0) (0) (1) (1) dent if we focus on the relative positioning of the pairs {T1 , T2 }, {T1 , T2 } (1) (1) and {T3 , T4 } with respect to a symmetry axis perpendicular to P placed at pz ≈ 2.7 (Fig. 1(a),2(a)). To the best of our knowledge, this is the first study to demonstrate the non-randomness of cumulative hydrophobic patterns around a NavCh’s pore. From an evolution point of view, the formation of membrane-parallel hydrophobic symmetries is of little surprise if we consider the channel and the membrane as coupled components of a larger system, i.e., the membrane-channel system, that has been optimized with respect to extra- and intracellular conditions. After all, pattern formation is a signaturephenomenon of non-equilibrium processes involving long-range spatial interactions [40]. Within this framework, the emergence and disappearance of hydrophobic pseudo-symmetries can be understood as an atomic self-organization phenomenon taking place alongside the structural transition from the PDs to the VSDs described by ν(p) (Fig. 1(a),2(a)). This implies that PDs and VSDs have evolved into distinct structural components not only for optimizing functional channel aspects such as voltage-sensing, but also for maximizing channel stability within the membrane. Given that the spatial transition from a proteins hydrophobic core to its hydrophilic exterior exhibits universal characteristics [37], our observations might be valid and generalizable to other NavChs, thus reflecting universal stability mechanisms. From a structural biology point of view, topological distributions φ(0) and (1) φ serve as molecular markers that can classify pore regions based on structural and hydrophobic channel properties. That was demonstrated was demonstrated by identifying locations of persisting extrema of φ(0) and φ(1) on P and reporting how these locations encode information about the functional modularity of the pore. In particular, the maxima of sm(φ(1) ) pinpoint locations on P which appear to correspond to the pore regions of the selectivity filter (SF) found at zSF ≈ −11.7, of the central cavity (CC) found at zCC ≈ 2.7 and of the AG found at zAG ≈ 22.7 (see also [19]). According to this scenario, pairing the maxima of sm(φ(0) ) with the minima of sm(φ(1) ) identifies locations on P where

11

a transition from one pore structural unit to the next takes place with respect ˆ: the transition from the extracellular funnel (EF) to the SF takes place for to z pz ∈ [−20.1, −17.9], from the SF to the CC for pz ∈ [−7.1, −5.3] and from the CC to the AG for pz ∈ [10.8, 18.0].

Figure 3: Summary of cumulative hydrophobic topological characteristics of the NavAb’s pore. In (a), (b), (c) and (d) we illustrate top-views (with respect to ˆ z) of the EF, of the SF, of the CC and of the AG, respectively. For clarity, only residue side chains sequences S180:V 184, T 175:W 179, 12

F 203:C217 and V 218:M 221 are illustrated, which are parts of the EF, of the SF, of the CC and of the AG, respectively. Atoms belonging to one of the porelining residue side chains M 181, S178, E177, L176, I210, V 213, C217, M 221 are colored according to a gradient that indicates the hydrophobic or the hydrophilic characteristics of the corresponding residue (SI Table S2, S3). Atoms that do not belong to any of the pore-lining residue side chains are shown in light grey color. (e), Illustration of the NavAb’s pore along P . For clarity, only atoms belonging to the PDs of subunit 1 and 3 are shown. (f), Traces of R(p), of sm(φ(0) ) and of sm(φ(1) ) for p ∈ Q. Pore points belonging to the subsets Q2·i+1,i=0,1,..,7 minimize their distance from the pore-lining residue side chains M 181, S178, E177, L176, I210, V 213, C217 and M 221, respectively (SI Table S3). Vertical (×) and (•) lines indicate locations on P where a minimum or maximum, respectively, of sm(φ(0) ) is found. Vertical (×) and (•) lines indicate locations on P where a minimum or maximum, respectively, of sm(φ(1) ) is found. ES stands for extracellular side, EF for extracellular funnel, SF for selectivity filter, CC for central cavity and AG for activation gate (see also [19]). Atomic illustrations were created with YASARA software (Yet Another Scientific Artificial Reality Application, www.yasara.org). Given that these locations account for a smoothed, multiscale hydrophobic effect, thermal fluctuations are expected to induce only small dislocations. In Fig. 3 we summarize these findings with respect to a simplified illustration of the NavAb’s pore. In the absence of radial structural channel symmetries a cumulative hydrophobic topological analysis is still possible by adding an additional analytical step, that is, the topological analysis of the spatial profile of the xy-component ~ (1) (p, lα (p)). For example, in the case of the heteromeric Nav1.7 channel, of m the xy-plane can be partitioned into four radial domains, each corresponding to a structural subunit, so that radial topological information can be extracted (1) by detecting the domain toward which the vector ~h xy (p, lα (p)) (SI equation S10) is pointing. Thus the presented algorithm can be immediately integrated in in-silico mutagenesis studies where a collection of human NavCh variants are selected for molecular screening (e.g., see [41]). In that case, a cumulative hydrophobic topological analysis might shed light upon atomic hydrophobic mechanisms that have a destabilizing effect on the current channel state. Eventually, this will improve already-existing tools for classification of NavChs variants of unknown pathophysiological significance and speed up the diagnostic pipeline. Moreover, it might open novel perspectives for molecular drugtherapeutic strategies via modulation of cumulative hydrophobic topological properties.

Methods 1. 3D structure preparation. PX4 ligands and water molecules were eliminated from the 3D structure of the NavAb (PDB code: 3RVY) and protonation 13

of the 3D structure was performed using the WHAT IF software [42, 43]. The principal axes of the protonated structure were estimated using the VMD software [44]. A global coordinate system (ˆ x, y ˆ, ˆ z) with origin O was introduced and the protonated structure was placed within it, so that the principal pore axis, i.e., the axis approximating the direction of the channel’s pore, is aligned with the z-axis. Orientation of the structure was set from the extracellular side (ES) to the intracellular side (IS) with respect to ˆ z and the channel’s molecular mass PNc 1 m ·c was set to coincide with O, where ci = (cx,i , cy,i , cz,i ) is center e = M i i i=1 the atomic center of the i-th atom, mi P is the mass of the i-th atom, Nc = 14776 Nc mi is the total molecular mass. is the total number of atoms and M = i=1 2. Geometrical characteristics of the pore. We consider P to represent the principal pore axis and p ∈ P to be a pore point. The pore radius at p, i.e., the radius of the smallest sphere that can be squeezed through the pore at p, is given by [45] (m1) R(p) = min {||ci − p|| − vdWi } i=1,2,..,Nc

, where ||·|| is the euclidean norm and vdWi is the Van der Waals radius of the i-th atom (SI Table S1). Consequently, the distance between p and its nearest neighbor atom is given by ¯ R(p) =

min

i=1,2,..,Nc

{||ci − p||}

(m2)

In analogy to equation m1, we introduce the outer surface radius at p L(p) =

max {||ci − p|| + vdWi }

i=1,2,..,Nc

(m3)

¯ The unit of measurement for R(p), R(p) and L(p) is expressed in [˚ A]. 3. Atomic sampling around P . In order to investigate cumulative atomic properties with respect to P , we introduced the p-dependent atomic sampling radius ( ¯ Nα ∈ Z+ L(p) − R(p) ¯ lα (p) = R(p) + α· for (m4) Nα 0 < α ≤ Nα , where Nα is the total number of nested sampling spheres centered at p and α denotes the index of the current sampling sphere. lα (p) plays the role of the molecular scale, i.e., it accounts for the size of the NavAb in ˚ A. The number of sampled atoms found within a sphere of radius lα (p) centered at p is given by n(p, lα (p)) =

Nc X

θ(lα (p) − ||ci − p||)

(m5)

i=1

, where θ(·) is the heaviside function. 4. Cumulative hydrophobic pore moments. Following [32] and [38], we introduced the zero-order cumulative hydrophobic pore moment h(0) (p, lα (p)) =

Nc X

θ(lα (p) − ||ci − p||)·HIξ,i

i=1

14

(m6)

, where HIξ,i = HIi + ξi is the i-th atomic hydrophobic index corresponding to the Kapcha-Rossky atomic hydrophobic scale [46] (SI Table S2) with additive noise ξi ∈ N (µ = 0, σ = 0.001) and the superscript (0) indicates the order of the moment and the measurement unit of h(0) (p, lα (p)) is roughly given by [kcal/mol] according to [32]. Moreover, we introduced the first-order, or dipole [32], cumulative hydrophobic pore moment ~h(1) (p, lα (p)) =

Nc X

θ(lα (p) − ||ci − p||)·HIξ,i ·~ ri,p

(m7)

i=1

~i,p is a vector from p to ci , the superscript (1) indicates the order , where r of the moment and the measurement unit of the magnitude ||~h(1) (p, lα (p))|| is roughly given by [kcal· ˚ A/mol]. In order to obtain an estimation of the bulk cumulative zero- and first-order hydrophobic pore moment effect originating from a collection of n(p, lα (p)) atoms, we introduced the zero-order mean cumulative hydrophobic pore moment m(0) (p, lα (p)) =

h(0) (p, lα (p)) n(p, lα (p))

(m8)

and the first-order mean cumulative hydrophobic pore moment ~ (1) (p, lα (p)) = m

~h(1) (p, lα (p)) n(p, lα (p))

(m9)

, respectively. 5. Discretization of equations (1)-(9). (a) Discretization along pdirection. We introduced the line grid Q = {p1 , p1 +∆p, .., pNp − ∆p, pNp } ⊂ P , where Np is the total number of grid pore points, ||∆p|| is the sampling distance between two consecutive pore points and p1 = (0, 0, pz,1 ), pNp = (0, 0, pz,Np ) are boundary pore points. Q is constructed by setting pz,1 = round( min (cz,i ), 1) = −27.1 and pz,Np = round( max (cz,i ), 1) = 26.8 i=1,2,..,Nc

i=1,2,..,Nc

with round(x ∈ R, n ∈ Z+ ) returning the value of x rounded upto the n-th decimal digit and by setting Np = 540 so that ||∆p|| = 0.1 ˚ A. The z-coordinate of an arbitrary p ∈ Q is represented as pz . (b) Discretization along lα (p)direction. We consider α ∈ A = {1, 2, .., Nα = 800}. 6. Detection of zero-crossing points of a scalar hydrophobic pore moment function along p-direction. We consider f (p, lα (p)) to represent a cumulative hydrophobic pore moment scalar function. If for a α ∈ A there is a pair {p0 = p − ∆p, p} ∈ Q for which the sign-change condition f (p0 , lα (p0 ))·f (p, lα (p)) < 0 is satisfied, then the 4D point (s = p0 +

|f (p0 , lα (p0 ))| ·∆p, lα (s)) + |f (p, lα (p))|

|f (p0 , lα (p0 ))|

(m10)

represents a zero-crossing point of f (p, lα (p)) along p-direction for a α ∈ A, 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 15

for a α ∈ A is represented as Γ(α). •

Author contributions M.N.X. conceived the presented idea and designed the algorithm; R.W., and P.L. helped with the refinement of algorithmic steps; M.N.X. performed the computations and data analysis; R.W., P.L., D.K., Y.Y., J.H., and S.G.W provided with critical feedback and helped with the interpretation of the results; Y.Y., and S.G.W encouraged M.N.X. to focus on specific aspects of the findings; B.H.S. helped with the supervision of the project; G.L. and C.G.F. were in charge of overall direction and planning; M.N.X. wrote the manuscript in consultation with R.W., P.L., J.H. and B.H.S.; PROPANE study group provided with technical and scientific support.

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

References [1] Hille, B. Ionic channels of excitable membranes, 3rd ed. (Sinauer Associates Inc., Sunderland MA, 2001). [2] Catterall, W.A. Voltage-gated sodium channels at 60: structure, function and pathophysiology. J. Physiol. 590, 2577-2589 (2012). [3] Cummins, T.R., Dib-Hajj, S.D., Waxman S.G. Electrophysiological properties of mutant Nav1.7 sodium channels in a painful inherited neuropathy. J. Neurosci. 24, 8232-8236 (2014). [4] Han, C. et al. Sporadic onset of erythermalgia: a gain-of-function mutation in Nav1.7. Ann. Neurol. 59, 553-558 (2006). [5] Harty, T.P. et al. N aV 1.7 mutant A863P in erythromelalgia: effects of altered activation and steady-state inactivation on excitability of nociceptive dorsal root ganglion neurons. J. Neurosci. 26, 12566-12575 (2006). [6] Lampert, A. et al. Erythromelalgia mutation L823R shifts activation and inactivation of threshold sodium channel Nav1.7 to hyperpolarized potentials. Biochem. Biophys. Res. Commun. 390, 319-324 (2006). [7] Stadler, T., O’Reilly, A.O., Lampert, A. Erythromelalgia mutation Q875E stabilizes the activated state of sodium channel Nav1.7. J. Biol. Chem. 290, 6316-6325 (2015). 16

[8] Yang, Y. et al. Mutations in SCN9A, encoding a sodium channel alpha subunit, in patients with primary erythermalgia. J. Med. Genet. 41, 171174 (2004). [9] Dib-Hajj, S.D. et al. Gain-of-function mutation in Nav1.7 in familial erythromelalgia induces bursting of sensory neurons. Brain. 128, 1847-1854 (2005). [10] Drenth, J.P. et al. SCN9A mutations define primary erythermalgia as a neuropathic disorder of voltage gated sodium channels. J. Invest. Dermatol. 124, 1333-1338 (2005). [11] Lee, M.J. et al. Characterization of a familial case with primary erythromelalgia from Taiwan. J. Neurol. 254, 210-214 (2007). [12] Drenth, J.P., Waxman, S.G. Mutations in sodium-channel gene SCN9a cause a spectrum of human genetic pain disorders. J. Clin. Invest. 117, 3603-3609 (2007). [13] Fertleman, C.R. et al. SCN9A mutations in paroxysmal extreme pain disorder: allelic variants underlie distinct channel defects and phenotypes. Neuron. 52, 767-774 (2006). [14] Jarecki, B.W., Sheets, P.L., Jackson, J.O. 2nd, Cummins, T.R. Paroxysmal extreme pain disorder mutations within the D3/S4-S5 linker of Nav1.7 cause moderate destabilization of fast inactivation. J. Physiol. 586, 41374153 (2008). [15] Dib-Hajj, S.D. et al. Paroxysmal extreme pain disorder M1627K mutation in human Nav1.7 renders DRG neurons hyperexcitable. Mol Pain. 4, 37 (2008). [16] Theile, J.W., Jarecki, B.W., Piekarz, A.D., Cummins, T.R. Nav1.7 mutations associated with paroxysmal extreme pain disorder, but not erythromelalgia, enhance Navβ4 peptide-mediated resurgent sodium currents. J. Physiol. 589, 597-608 (2011). [17] Faber, C.G. et al. Gain of function Nav1.7 mutations in idiopathic small fiber neuropathy. Ann. Neurol. 71, 26-39 (2012). [18] Hoeijmakers, J.G.J. et al. Small nerve fibres, small hands and small feet: a new syndrome of pain, dysautonomia and acromesomelia in a kindred with a novel NaV1.7 mutation. Brain. 135, 345-358 (2012). [19] Payandeh, J., Scheuer, T., Zheng, N., & Catterall, W.A. The crystal structure of a voltage-gated sodium channel. Nat. 475, 353-358 (2011). [20] Yonkunas, M., & Kurnikova, M. The hydrophobic effect contributes to the closed state of a simplified ion channel through a conserved Hydrophobic patch at the pore-helix crossing. Front. Pharmacol. 6, 284 (2015). 17

[21] Kitaguchi, T., Sukhareva, M., & Swartz, K.J. Stabilizing the closed S6 gate in the Shaker Kv channel through modification of a hydrophobic seal. J. Gen. Physiol. 124, 319-332 (2004). [22] Beckstein, O., Sansom, M.S. Liquid-vapor oscillations of water in hydrophobic nanopores. Proc. Natl. Acad. Sci. U.S.A. 100, 7063-7068 (2003). [23] Beckstein, O., & Sansom, M.S. The influence of geometry, surface character, and flexibility on the permeation of ions and water through biological pores. Phys. Biol. 1, 42-52 (2004). [24] Allen, R., & Hansen, J.-P. Molecular dynamics investigation of water permeation through nanopores. J. Chem. Phys. 119, 3905-3919 (2003). [25] Jensen, M.Ø. et al. Principles of conduction and hydrophobic gating in K + channels. Proc. Natl. Acad. Sci. U.S.A. 107, 5833-5838 (2010). [26] Lampert, A. et al. A pore-blocking hydrophobic motif at the cytoplasmic aperture of the closed-state Nav1.7 channel is disrupted by the erythromelalgia-associated F1449V mutation. J. Biol. Chem. 283, 2411824127 (2008). [27] Yang, Y., Estacion, M., Dib-Hajj, S.D., & Waxman S.G. Molecular architecture of a sodium channel S6 Helix: radial tuning of the voltage-gated sodium channel 1.7 activation gate. J. Biol. Chem. 288, 13741-13747 (2013). [28] Hille, B. Local anesthetics: hydrophilic and hydrophobic pathways for the drug-receptor reaction. J. Gen. Physiol. 69, 497-515 (1977). [29] Israelachvili, J., & Pashley, R.M. The hydrophobic interaction is long range, decaying exponentially with distance. Nat. 300, 341-342 (1982). [30] Hammer, M.U., Anderson, T.H., Chaimovich, A., Shell, M.S., & Israelachvili, J. The search for the hydrophobic force law. Faraday Discuss. 146, 299-308 (2010). [31] Tabor, R.F., Grieser, F., Dagastine, R.R., & Chan D.Y.C. The hydrophobic force: measurements and methods. Phys. Chem. Chem. Phys. 16, 1806518075 (2014). [32] Eisenberg, D., Weiss, R.M., Terwilliger, T.C., & Wilcox, W. Hydrophobic moments and protein structure. Faraday Symp. Chem. Soc. 17, 109-120 (1982). [33] Reißer, S., Strandberg, E., Steinbrecher, T., & Ulrich A.S. 3D hydrophobic moment vectors as a tool to characterize the surface polarity of amphiphilic peptides. Biophys. J. 106, 2385-2394 (2014). [34] Doyle, D.A. et al. The structure of the potassium channel: Molecular basis of K + conduction and selectivity. Science 280, 69-77 (1998).

18

[35] Meirovitch, H., & Scheraga, H.A. Empirical studies of hydrophobicity. 3. Radial distribution of clusters of hydrophobic and hydrophilic aminoacids. Macromolecules 14, 340-345 (1981). [36] Silverman, B.D. Hydrophobic moments of protein structures: Spatially profiling the distribution. Proc. Natl. Acad. Sci. U.S.A. 98, 4996-5001 (2001). [37] Silverman, B.D. Spatial profiling of protein hydrophobicity: Native vs. decoy structures. Proteins 52, 561-572 (2003). [38] Silverman, B.D. Hydrophobicity of transmembrane proteins: Spatially profiling the distribution. Protein Sci. 12, 586-599 (2003). [39] Silverman, B.D. Three-dimensional moments of molecular property fields. J. Chem. Inf. Comput. Sci. 40, 1470-1476 (2000). [40] Gollub, J.P. & Langer, J.S. Pattern formation in nonequilibrium physics. Rev. Mod. Phys. 71, S396 (1999). [41] Kapetis, D. et al. Network topology of NaV1.7 mutations in sodium channel-related painful disorders. BMC Systems Biology 11, 28 (2017). [42] Vriend, G. WHAT IF: A molecular modeling and drug design program. J. Mol. Graph. 8, 52-56 (1980). [43] Hooft, R.W., Sander, C., Vriend, G. Positioning hydrogen atoms by optimizing hydrogen-bond networks in protein structures. PROTEINS 26, 363376 (1996). [44] Humphrey, W., Dalke, A., & Schulten, K. VMD-Visual Molecular Dynamics. J. Molec. Graph. 14, 33-38 (1996). [45] Smart, O.S., Neduvelil, J.G., Wang, X., Wallace, B.A., & Sansom M.S. HOLE: A program for the analysis of the pore dimensions of ion channel structural models. J. Mol. Graph. 14, 354-360 (1996). [46] Kapcha, L.H., & Rossky, P.J. A simple atomic-level hydrophobicity scale reveals protein interfacial structure. J. Mol. Biol. 426, 484-498 (2014). [47] R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria (2014) http://www. R-project.org/.

19