Hydropathicity-based prediction of pain-causing NaV1.7 variants.

BACKGROUND: Mutation-induced variations in the functional architecture of the NaV1.7 channel protein are causally related to a broad spectrum of human pain disorders. Predicting in silico the phenotype of NaV1.7 variant is of major clinical importance; it can aid in reducing costs of in vitro pathophysiological characterization of NaV1.7 variants, as well as, in the design of drug agents for counteracting pain-disease symptoms.
RESULTS: In this work, we utilize spatial complexity of hydropathic effects toward predicting which NaV1.7 variants cause pain (and which are neutral) based on the location of corresponding mutation sites within the NaV1.7 structure. For that, we analyze topological and scaling hydropathic characteristics of the atomic environment around NaV1.7's pore and probe their spatial correlation with mutation sites. We show that pain-related mutation sites occupy structural locations in proximity to a hydrophobic patch lining the pore while clustering at a critical hydropathic-interactions distance from the selectivity filter (SF). Taken together, these observations can differentiate pain-related NaV1.7 variants from neutral ones, i.e., NaV1.7 variants not causing pain disease, with 80.5[Formula: see text] sensitivity and 93.7[Formula: see text] specificity [area under the receiver operating characteristics curve = 0.872].
CONCLUSIONS: Our findings suggest that maintaining hydrophobic NaV1.7 interior intact, as well as, a finely-tuned (dictated by hydropathic interactions) distance from the SF might be necessary molecular conditions for physiological NaV1.7 functioning. The main advantage for using the presented predictive scheme is its negligible computational cost, as well as, hydropathicity-based biophysical rationalization.

Introduction

Voltage-gated sodium channels (NaVChs) are pore-forming proteins embedded in cell membranes. They are members of the ion channels superfamily and their main physiological role is to control transport of sodium ions into the cell. The human NaV1.7 channel is encoded by the SCN9A gene and is preferentially expressed in peripheral neurons (e.g., dorsal root ganglion (DRG) nociceptors) responsible for networking pain signals. The structure of the NaV1.7 -subunits is that of a pore-forming tetramer via assembly of four heterogeneous domains (DI-DIV). Three intracellular linkers (L1-L3) form structural interconnections among subsequent domains. Each domain comprises six transmembrane helices (S1-S6) organized into a pore module (PM) forming an ion-conduction pathway coupled with a voltage-senor (VS). Mechanistic description of NaV1.7’s function is that VSs react to extracellular changes in ionic concentrations by moving outwards thus exerting a pulling force upon the PM which opens the channel pore. Closed-to-open conformational change leads to channel activation, i.e., renders it conductive to sodium ions. Missense mutations in the SCN9A gene can destabilize the NaV1.7’s functional architecture thus disrupting physiological sodium-ions current and, consequently, deregulate flow of sodium ions through the pore. At a cellular level, these genetically-caused destabilizations can affect neuronal excitability by inducing a gain-of-function (GOF) effect, i.e., by increasing the net sodium-ions currents, thus triggering a wide spectrum of pain diseases such as inherited erythromelalgia (IEM) [1-30], paroxysmal extreme pain disorder (PEPD) [31-37] and small fiber neuropathy (SFN) [38-42]. A proof of concept for the GOF-pain correlation hypothesis came from identification of missense SCN9A-gene mutations inducing a loss-of-function (LOF) effect, i.e., decreasing net sodium-ions current, that is causally related to clinical symptoms of loss of pain sensation [43-45].

Hydropathic interactions (HIs) represent a summary of fundamental molecular interactions [46] driving molecular phenomena such as protein folding, protein hydrophobic-core stability, self-assembly of amphiphilic molecules, and “dewetting” of molecular surfaces (for a review in HIs-driven phenomena see [47]). Within the field of ion channels research, experimental and computational studies have shown that pore-lining patches of hydrophobic residues are crucial regulators of, both, pore’s gating behavior (a phenomenon termed as “hydrophobic gating” [48]) and channel stability [49], e.g., via formation of hydrogen-bonds networks expanding through their surroundings [50, 51]. Crucially, hydrophobic patches (HPs) are widely conserved across voltage-gated channels and often associated with the formation of centrally-located cavity (CC) [52]. Mutations perturbing this hydrophobic motif can lead to drastic changes of net ion currents [53-56].

Computational modeling of HIs combined with biophysical observations extracted from in vitro NaVCh pathophysiological characterization can propel our understanding of mechanistic linkages between mutation-induced perturbations and human pain pathophysiology. Key-studies toward this direction were these of Lampert et al [57], and of Yang et al [58] demonstrating how the F1449V mutation, and the in-frame-deletion L955Del, respectively, can disrupt a hydrophobic ring stabilizing the putative activation gate (AG) of the NaV1.7 thus acting as disease-causing molecular triggers. Moreover, computational modeling successfully deduced an energetic coupling between two different IEM-related mutations foreseen by their geometrical proximity in NaV1.7 structure [7]. A question that naturally arose from these studies was whether a detailed examination of HIs network characteristics within a NaV1.7 structure can reveal statistically-significant but also biophysically-relevant differentiations among the WT structure and its variants. This question was probed by Kapetis et al [59]; a network-theoretical computational framework was introduced in order to capture changes in interatomic HIs within a NaV1.7 structural model induced by pain-related mutations. The study reported on a betweenness centrality network measure achieving a statistically-important differentiation of pain-related variants from a collection of neutrals, i.e., variants not causing pain disease. Notably, this approach highlighted the prominent role that HIs play in NaV1.7’s stability and reported on plausible mutation-mechanism scenarios disrupting hydrophobic contacts among neighboring and distant residues. Another remark on the multi-scale nature of HIs was made by the authors of [60] suggesting that a pathogenic mutation in the KCNA1 gene encoding the human voltage-gated potassium channel KV1.1 can de-tune HIs equilibrium (and, consequently, destabilize KV1.1’s pre-open conformation) implying that mutation-induced perturbation effects can destroy finely-tuned network-like HIs expanding throughout the structure as a whole. Interestingly, the fine-tuning hypothesis was proposed also for the NaV1.7; a recent study employing a machine learning (MLE) computational pipeline for predicting NaV1.7’s variant pathogenicity suggested that the fine-tuning of the NaV1.7 channel is so delicate that limits classification accuracy of practically any computational approach [61]. Taken together, these observations highlight the highly-cooperative nature of HIs [46] and suggest that even small changes in the hydropathic spatial distribution profile of a channel structure can have a detrimental impact on the functional architecture which, in turn, might induce clinically-observed alterations of electrophysiology.

Following [59, 61], this study aims at probing the finely-tuned hypothesis for the NaV1.7’s atomic hydropathic environment towards predicting whether a NaV1.7 variant causes pain or not. We demonstrate our approach on a closed-state NaV1.7 structural model (first presented in [7] and later also used in [30, 58]) by investigating cumulative, i.e., scale-dependent, hydropathic properties of its porous atomic environment in relation to structural locations of missense SCN9A-gene mutations. In order to tackle spatial complexities emerging from the highly-cooperative nature of HIs we adopt a modeling approach rooted in the hypothesis that proteins can be represented as self-organized criticality (SOC) [62] archetypes; protein structures are thought to have been evolutionary optimized with respect to extrema in some thermodynamic property (or properties) capturing a qualitative reorganization of the atomic environment [63, 64]. The intra-channel locations where these macroscopic thermodynamic changes take place correspond to so-called critical points of the atomic structure [63, 64]. The highly-cooperative nature of HIs has placed structure-retrieved hydropathic properties in the epicenter of SOC hypothesis [63-66]). It is important to note that computational evidence for a universal hydrophobic-to-hydrophilic (or inside-outside with “inside” referring to the hydrophobic core, and “outside” referring to the hydrophilic exterior) spatial transition in protein systems was first provided before the formulation of the SOC hypothesis (see [67]). Departing from this phenomenological basis, we here utilize the finite-size scaling analysis methodologies presented in [68, 69] for screening hydropathic morphology around NaV1.7’s pore [68] toward identification of critical points associated with NaV1.7’s functional architecture. Biophysical relevance of retrieved observations is justified not only in terms of the scale-invariance of a carefully-chosen cumulative hydropathicity-property function but also with respect to conserved structural NaV features such as the PM-VSs spatial transition [70, 71] and the architecture of the selectivity filter (SF) [72]. In particular, we demonstrate that the atomic cumulative distribution function around NaV1.7’s pore exhibits a sigmoid profile with inflection points matching closely the conserved PM-VSs spatial transition. This provides a rigorous description of atom-packing geometry and, consequently, a macroscopic partitioning of the atomic environment around the pore allowing for mapping the spatial profile of the atomic cumulative hydropathicity-property function and mutation sites on two dimensions. The SOC hypothesis is then accepted (or rejected) depending on whether the cumulative hydropathicity-property function is globally maximized and exhibits power-law-like scaling behavior in the vicinity of the inflection point (or not).

Our mapping procedures reveal the formation of a HP incorporating NaV1.7’s CC and AG. We report on two “hot” map areas attracting pain-related mutation sites distributed along HP’s boundary. Probing the SOC hypothesis for HIs around NaV1.7’s pore reveals that “hot” structural locations tend to cluster at a distance of 33.4 Å from the SF. Stability implications of these observations can be deduced by considering that in the vicinity of the critical point the range and intensity of HIs increase in a power-law fashion thus favoring amplification and propagation of mutation-induced perturbations peripherally to the HP and at critical HIs-distance from the SF thus not directly affecting neither of them. The clinical translational value of our findings is tested by predicting pathogenicity of 84 NaV1.7 variants; a weighted average of HP- and SF-related distance metrics can classify up to 29 (out of 36) pain-related variants and 45 (out of 48) neutral variants correctly.

Methods

All computations were performed in R [73] environment unless stated differently.

3D structure preparation

The NaV1.7 atomic structural model used throughout this study was constructed via homology modeling procedures based on the pre-open NaVAb template [PDB code: 3RVY] [74] according to [7]. In short, the first step undertaken was to generate structural models of the four transmembrane domains (DI-DIV) by utilizing the membrane-bound protein predication algorithm GPCR-ITASSER [75-77]. Then, each domain was aligned upon the pre-open NaVAb template by using the TM-align algorithm [78]. Finally, the four transmembrane domains were placed together in a clockwise order viewed from extracellular side (ES) according to [79, 80], and the resulting four domain complex structural model was refined by fragment-guided molecular dynamics simulations aiming at removing interdomain clashes, optimizing hydrogen atoms network, and, consequently, increasing model quality [75, 76, 81]. The retrieved model consists of 1140 protonated amino acids (DI: P229:K417, DII: P839:T972, DIII: M1296:G1461, DIV: K1617:T1763). A comparison via superposition of the retrieved model with a recently resolved cryo-electron microscopy (cryo-EM) NaV1.7 structure [PDB code: 6J8J] [82] is presented in Additional file 1: S1.

In continuation, principal axes of the NaV1.7 model structure were estimated by using the VMD software [83]. A global coordinate system was introduced with its center at O and the NaV1.7’s principal pore axis, i.e., the axis approximating the direction of the channel’s pore, was aligned with the z-axis with orientation from the ES toward the intracellular side (IS) with respect to . The atomic center of the 3D structure was set to overlap with O, where is the atomic center of the ith atom, is the mass of the ith atom, is the total number of atoms, and is the total molecular mass (values of atomic masses are the same as [68, 69]).

Geometrical characteristics of the pore

We navigated through the skewed NaV1.7’s pore by introducing a set of pore points (see Additional file 1: S2). The pore radius at is given by [? ]where is the euclidean norm and vdW is the van der Waals radius of the ith atom (values of van der Waals radii are the same as [68, 69]). The distance between and its nearest neighbor atom corresponds then toand the outer surface radius at is given by [68, 69]where the unit of measurement for , and is expressed in [Å].

Finite-size sampling around the pore

The atomic environment around is sampled with concentric spheres placed at of increasing radius [68, 69]where is the total number of sampling spheres and denotes the index of the sampling sphere (i.e., the scaling index). indicates the size, i.e., molecular scale, of the spherical cluster of atoms around in [Å] and the finite channel size measured with respect to . Accordingly, the atomic cumulative distribution function (CDF) at is given by [68, 69]where is the heaviside function. Note that essentially describes how atoms are packed around . In computational practice is set to be “large enough” approximating the continuous case via dense sampling.

Mathematical modeling of atomic CDF

Modeling of the atomic CDF was performed by employing the GROFIT routine [84]. Hypothesis underlying this modeling approach is that the atomic CDF for a given follows a sigmoid pattern. Hypothesis is approved if GROFIT manages to fit any of available candidate models to the CDF data. Available candidate models are a re-parametrized algebraic forms [85] of the Logistic model [86], of the Gompertz model [87]with , of the the modified Gompertz model [85]and of the Richards model [88]were fitted on along -direction where are model parameters. The candidate model that best fits the atomic CDF is then selected based on minimization of an Akaike information criterion (see [84] for algorithmic details).

Following [69], model parameters interpretation was performed with respect to the inflection pointthat determines the location along -direction where the radial distribution function (RDF), , maximizes. The RDF maximum value is given by the parameter accounting for the maximum atom-packing rate (or, equivalently, for the maximum atomic density) around . Parameter is the asymptote value of the fitted model, i.e., , describing what happens when becomes arbitrary large. Parameter determines the location along -direction where the lag atom-packing domain ends, i.e., its size. Notably, parameter can be expressed in terms of the ratio with determining the location along -direction where the asymptote atom-packing domain begins. Parameter affects the shape of the Richards model curve, as well as, the location of the inflection point thus plays the role of the summary atom-packing parameter. Parameters , and of the modified Gompertz model indicate the location, and the slope, respectively, of a second increase in the modified Gompertz model curve [84]. The Logistic and the Gompertz model are retrieved from the Richards model for , and for , respectively, as shown in [89], thus they are considered as special cases of the Richards model.

Cumulative hydropathicity-property functions

The hydropathic density of the atomic environment around was approximated in terms of [68]where corresponds to the cumulative zero-order hydropathic pore moment function [90] with representing the ith atomic hydrophobic index in accordance with the Kapcha&Rossky atomic hydropathic scale [91] with additive Gaussian noise . The superscript “(0)” indicates the moment order.

The hydropathic interatomic interaction strength (HIIS) at , i.e., the average interaction strength between an atomic component found within the cluster of size and its surroundings, was approximated in terms of the hydropathic imbalance (or interaction strength) pore function [68, 69]where corresponds to the cumulative first-order hydropathic pore moment function [90] quantifying the hydropathic inter-cluster interaction strength (HIcIS) with being the vector from to . The superscript “(1)” indicates the moment order.

Introduction of a weak noise source practically guarantees that scalars and are non-zero for every combination of and while their scaling behavior remains practically unaffected. Throughout this study we consider pore’s physichochemical field characteristics to be adequately described in terms of the HIIS axial field component, given that the magnitude of the radial field component, , is statistically expected to remain always smaller than after a cut-off, lag-domain scale and, thereafter, to gradually shrink (see Additional file 1: S3). Accordingly, we focus only on the scaling behavior and topology of the HIIS axial field component which can occupy only two orientation-states; an “in” orientation-state which is characterized by pointing towards the intracellular side (IS), i.e., by , and an “out” orientation-state which is characterized by pointing towards the ES, i.e., by . Topological changes in HIIS axial field component are detected according to the algorithmic scheme presented in [68] (see Additional file 1: S4).

Finite-size scaling of HIIS axial field component

In accordance to [69], a scale-invariant interval of the HIIS axial field component corresponds to combinations of with for which the power-law approximationis accurately satisfied indicating that HIIS axial field component stabilizing the cluster of atoms around span a range up to Å. Sign of quantifies the rate at which intensity and range of HIIS axial field component increase or decrease for increasing atomic cluster size. From a HIs-network standpoint, indicates whether HIs network interconnectivity is up- or down-regulated, i.e., whether HIs contacts, e.g., hydrogen bonds, are created or destroyed within the structure. The energy levels associated with HIs contacts creation (or destruction) guaranteeing stability of the atomic cluster can then be approximated bymeasured in kcal/(molatomic cluster).

Results

In Fig. 1 we demonstrate that atomic CDF around NaV1.7’s pore follows a sigmoid profile which can be adequately described by the Richards model (Additional file 1: S5). The atomic environment around the pore can thus be partitioned into three consecutive atom-packing domains spanning the channel from the inside (i.e., pore-lining environment) to the outside (i.e., voltage-sensing environment) (Fig. 1b, c). The innermost domain corresponds to the lag domain consisting primarily of structural elements lining the pore (approximately 95 of lag-domain structural components belong to S5–S6 helices). The outermost domain corresponds to the asymptote domain consisting mostly of structural elements drawn from the S1–S4 voltage-sensing helices (approximately 63.3 of asymptote-domain structural components belong to the S1–S4 helices). In-between the lag and asymptote domain, a structurally-heterogeneous inflection domain is formed consisting of two parts separated by inflection points, . Structural locations of the inflection points correspond to intra-channel regions where the atomic RDF maximizes and, as demonstrated in Fig. 1a, they closely match the PMs-VSs spatial transition described by (see Additional file 1: S6 for calculation of PMs-VSs spatial transition characteristics). Accordingly, serves here as a macroscopic boundary line splitting the atomic environment around NaV1.7’s pore into two phases, namely, a pre-inflection phase for and a post-inflection phase for accounting for atomic sub-environments containing mainly structural components belonging to the PMs and VSs, respectively (Fig. 1b).

Fig. 1

Atom-packing around NaV1.7’s pore. a Cartoon illustration of the NaV1.7 structural model (side view). b Cartoon illustration of the NaV1.7 structural model (intracellular-to-extracellular view). The atomic environment around every pore point is partitioned into three consecutive atom-packing domains; a lag domain realized for , an inflection domain consisting of two parts realized for and , respectively, and an asymptote domain realized for (see “Methods” section). , and are represented in b in terms of , , and , respectively, roughly indicating the median-statistical value of the radial distance from at which the transition among subsequent domains takes place. The median-statistical value of the radial distance from at which the PMs-to-VSs transition, , takes place is also illustrated. Note that and are almost indistinguishable, i.e., . c Traces of statistical descriptors of the normalized (with respect to ) atomic CDF, , and of its best-fitted Richards model are plotted in log-scale with shaded areas around indicative of 95 confidence intervals. , , and are represented in c in terms of statistical descriptors , , and , respectively, returning the median-statistical value of corresponding scaling indices . All statistical descriptors correspond to median values and are calculated according to the scheme presented in Additional file 1: S7

Based on the geometrical partition scheme summarized in Fig. 1 we proceeded with mapping of missense SCN9A-gene mutations on two-dimensions based on their intra-channel structural locations. For that, we utilized a collection of well-studied GOF NaV1.7 mutations phenotypically related with IEM, SFN and PEPD pain disease (total number of pain-related mutation sites: 36) (Additional file 1: S8). In addition, a collection of neutrals SCN9A-gene mutations was introduced consisting of SCN9A-gene homologous single-nucleotide polymorphisms (hSNPs) and of SCN9A-gene variants not causing biophysical abnormalities (nBABNs) (imported from [59]), as well as, of a small number of SCN9A-gene variants previously classified as non-pathogenic based on algorithmic procedures described in [92] (total number of neutral mutation sites: 48) (Additional file 1: S8). Given that neutrals are not expected to associate with pain disease phenotypes, they are less likely to trigger any significant NaV1.7 destabilizations and/or to perturb conserved hydropathic characteristics around NaV1.7’s pore.

The majority (i.e., ) of pain-related mutation sites are distributed within the inflection domain with a tendency to cluster toward a centrally-located map area along the inflection-points line (see area “a2” on Fig. 2a in relation to Fig. 2b, c). On the other hand, the majority (i.e., ) of neutral mutation sites are distributed within the second part of the inflection domain so that their map occupancy rate tends to maximize approximately in the middle of the second part of the inflection domain (see area “a3” on Fig. 2a in relation to Fig. 2b, c). The rest of pain-related mutation sites are found within the lag domain toward the intracellular side of the NaV1.7. In particular, their map occupancy rate maximizes in the vicinity of the boundary line toward the IS (see area “a1” on Fig. 2a in relation to Fig. 2b, c). Taken together, these observations indicate that for decreasing molecular scale the probability of a missense SCN9A-gene mutation to translate into a pain-related phenotype increases which is of little surprise considering that mutations perturbing NaV1.7’s interior are more likely to perturb tight packing of S5-S6 helices and, consequently, produce electrophysiological alternations.

Fig. 2

Mapping of missense SCN9A-gene mutation sites around NaV1.7’s pore. a Mapping of missense SCN9A-gene mutation sites for and . Two sets of missense SCN9A-gene mutation sites are employed; a pain-related set containing IEM, PPD and SFN mutation sites, and a neutral set containing mutation sites which are not expected to associate with pain disease phenotypes (Additional file 1: S8). Scaling indices lines , and highlight the boundaries among consecutive atom-packing domains. Specifically, denotes the ending and beginning of the lag and inflection domain, respectively. denotes the ending and beginning of the inflection and asymptote domain, respectively. denotes the location of inflection points and the ending and beginning of the first and second part of the inflection domain, respectively. Labels “a1”, “a2”, and “a3” indicate map areas where the number of mutation sites maximizes, i.e., mutation sites occupancy rates maximize. a Map occupancy rate of mutation sites along -direction. c Map occupancy rate of mutation sites along -direction. Red- and blue-colored histograms account for map occupancy rates of pain-related, and neutral mutation sites, respectively

Mapping of hydropathic density profile reveals the formation of a HP lining NaV1.7’s pore (see Fig. 3, blue-colored domains and ). Specifically, we demonstrate that the center of the pore is lined by predominantly hydrophobic atomic components expanding toward the IS where occlusion of the pore takes place by the ring of Y405 (DI), F960 (DII) F1449 (DIII) and F1752 (DIV) residues which are known to form the NaV1.7’s activation gate (AG) (see Fig. 3, blue-colored domain ). Hydrophobic expansion toward the ES is disrupted by the weakly-hydrophilic SF so that the pore is lined by a small-sized hydrophobic cluster located in-between the strongly-hydrophilic ES mouth and the weakly-hydrophilic SF (see Fig. 3, blue-colored domain ). Macroscopically, the wide CC translates into a structural contraction event as the outer surface radius is locally minimized so that the channel is split into two funnel-like structural compartments (see Fig. 3, trace of ). Structural combination of the PMs with the VSs results in a hydrophilic atomic environment as indicated by the red-colored domain covering the largest map area and incorporating both the SF and the ES mouth, as well as, the IS channel end (see Fig. 3). The SF’s microenvironment is formed by the residues D361 (DI), E930 (DII), K1406 (DIII) and A1698 (DIV) where a bare sodium ion of radius Å  can exactly fit in (Fig. 3).

Fig. 3

Spatial profile of the hydropathic density around NaV1.7’s pore. Contour map of the hydropathic density pore function, , is illustrated for and . Blue- and red-colored contour domains represent hydrophobic and hydrophilic domains around the pore, respectively. Black lines , and indicate geometrical pore characteristics (see “Methods” section). Magenta dashed line depicts the scales at which the PMs-VSs spatial transition takes place. Dashed black lines , and account for the boundaries among subsequent atom-packing domains (see “Methods” section). Zero-crossing points of collected in define HP’s boundary, i.e., the boundary between HP-forming domains and , and hydrophilic domain (see Additional file 1: S4 for calculation of zero-crossing points and construction of ). Black arrows , , and highlight domain boundaries. Two sets of missense SCN9A-gene mutation sites are employed; a pain-related set containing IEM, PPD and SFN mutation sites, and a neutral set containing mutation sites which are not expected to associate with pain disease phenotypes (Additional file 1: S8). Mutation sites highlighted in red color correspond to misclassified events (classification criterion; median distance from HP’s boundary (see Additional file 1: S9a)). Grey-shaded areas “a1”, “a2”, and “a3” highlight contour map regions where the number of mutation sites maximizes, i.e., mutation sites occupancy rates maximize. ES, SF, CC, AG, and IS labels mark the locations of the extracellular side, of the selectivity filter, of the central cavity, of the activation gate, and of the intracellular side, respectively

Strikingly, approximately of pain-related mutation sites are found within with a tendency to cluster along HP’s boundary as map areas “a1” and “a2” occupy contour map regions where the transition from to takes place (Fig. 3). On the other hand, of neutral mutation sites are located within the and only one neutral mutation site found within , while area “a3” is distributed solely within (Fig. 3).

Given that mutations perturbing an ion channel’s hydrophobic interior properties pose a high risk for ion current alternations [53-56], we hypothesize that mutations occurring at structural locations in proximity to the HP are more likely to be related with GOF pain phenotypes. We tested this hypothesis by statistically approximating the distance between each mutation structural location and HP’s boundary (Additional file 1: S9a) and feeding retrieved median distances into a binary classifier. We achieved to classify correctly 29 (out of 36) and 38 (out of 48) of pain-related and neutral, respectively, mutations correctly with a cut-off median distance of 18.13 Å  (Fig. 4). This translates to an area under receiver operating characteristics (ROC) curve of 0.787 and pain phenotype prediction with specificity of 0.791 and sensitivity of 0.805 (Fig. 4a).

Fig. 4

Binary classification of missense SCN9A-gene mutation sites based on their median distance from HP’s boundary. a ROC-curve plot constructed from data of median distances between mutation sites and HP’s boundary (for construction of data set see Additional file 1: S9a). Optimal threshold value corresponds to specificity and sensitivity values of 0.791 and 0.805, respectively. Area under ROC curve is 0.787. b Visualization of ROC curve data. Optimal threshold value 18.13 Å is marked with black dashed line. Shaded area around median distance values indicates the 95 confidence intervals. ROC curve is constructed in R [73] by using the pROC package [93]. Two sets of missense SCN9A-gene mutation sites are employed; a pain-related set containing IEM, PPD and SFN mutation sites, and a neutral set containing mutation sites which are not expected to associate with pain disease phenotypes (Additional file 1: S8)

Misclassified pain-related mutations S211P, L823R, W1538R, I720K, I739V and T1596I are found outside of the HP thus not in hydrophobic proximity neither to CC nor to AG (Figs. 3 and 4b). Only a single pain-related misclassification is found within HP, namely, R185H that is located with (Figs. 3 and 4b). This striking misclassification is due to the elementary statistical approach adopted for calculating distance scores which fails to fully capture the complex topology of HP-forming and (Additional file 1: S8). Misclassified neutral mutations V1428I, T920N, V194I, V1613I, T1398N, I1399D, S1419N, D1662A, D1674A and K1700A are found inside the HP, and, more precisely, in proximity to HP’s boundary, with a tendency to cluster around the SF (Figs. 3 and 4b).

We further investigated the relation between cumulative hydropathic characteristics and mutation sites around NaV1.7’s pore by focusing on the HIIS axial field component illustrated in Fig. 5a. As we demonstrate in in Fig. 5a, HIIS axial field component topology can be summarized into five domains, namely, , , , and (Fig. 5a). The centrally-located domain covers the largest map area and roughly dichotomizes the contour map into two pseudo-symmetric parts, namely, an ES part incorporating and , and an IS part incorporating and . Pain-related mutation sites are solely found within the and domains with occupancy rates of 58 and 42, respectively. On the other hand, neutral sites are found within the , and with occupancy rates of 14, 19 and 67, respectively.

Fig. 5

Spatial profile of HIIS axial field component along NaV1.7’s pore. a Contour map of HIIS axial field component, , for and . Blue- and red-colored contour domains represent configurations of with orientation “out” and “in”, respectively. Black lines , and indicate geometrical pore characteristics (see “Methods” section). Magenta dashed line depicts the scales at which the PMs-VSs spatial transition takes place. Dashed black lines and account for the upper and lower boundary of the lag and asymptote atom-packing domains. Critical radius Å is plotted while inflection points line is omitted for clarity. Zero-crossing points of collected in (see Additional file 1: S4) correspond to boundaries among contour domains , , , and . Black arrows , , , , and are plotted in order to highlight domain boundaries. Two sets of missense SCN9A-gene mutation sites are employed; a pain-related set containing IEM, PPD and SFN mutation sites, and a neutral set containing mutation sites which are not expected to associate with pain disease phenotypes (Additional file 1: S8). Mutation sites highlighted with red color correspond to misclassified events (classification criterion; distance from the SF (see Additional file 1: S9b)). Grey-shaded areas “a1”, “a2”, and “a3” highlight contour map regions where the number of mutation sites maximizes, i.e., mutation sites occupancy rates maximize. ES, SF, CC, AG, and IS labels mark the locations of the extracellular side, of the critical pore point , of the central cavity, of the activation gate, and of the intracellular side, respectively. b Traces of normalized (with respect to corresponding maximum values) median distances between pain-related and neutral mutation sites from the critical point, , represented as and , respectively, are plotted for (see Additional file 1: S9b for calculation of and ). Circles indicate that for , is globally minimized with , while exhibits a local maximum with . c Trace of for and . Power-law approximations of described by Eq. r1 are plotted in light and dark green color accounting for the first and second part of the inflection domain, i.e., for and , respectively. The mean absolute relative fitting errors (MARFEs) of the power-law approximation for the first and second part of the inflection domain are and , respectively. d Trace of AE, , for and . Modeling approximations of described by Eq. r2 are plotted in light and dark green color accounting for the first and second part of the inflection domain, i.e., for and , respectively. The MARFEs of the modeling approximation for the first and second part of the inflection domain are and , respectively. Extrapolation of model approximations toward the lag domain, i.e., for , and toward the asymptote domain, i.e., for , are plotted with dashed light and dark green lines, respectively, and result in a MARFE of and , respectively. Richards model parameters used for modeling AE are

In order to decode mutation sites clustering behavior on the contour map of Fig. 5a we adopted a phenomenological approach that presumes the existence of a critical point, , associated with the spatial organization of HIs around the SF [63, 64, 69] with being a critical SF pore point coordinate (Additional file 1: S10). A crucial result that motivated us to adopt such an approach is that pain-related mutation sites are attracted toward the critical point in sheer contrast to neutral mutation sites which are repelled from it (Fig. 5b). We refer to this phenomenon with the term critical clustering. Geometrically, the formation of the critical mutation-sites cluster reflects the tendency of structural locations of pain-related mutations to minimize their distance from the surface of the sphere of radius Å  around the SF; intuitive graphical representation of this phenomenon is provided in Fig. 5a where we show that “hot” areas I and II intersect with the green radius line representing critical sphere’s surface.

The scaling behavior of HIIS axial field component around is adequately described in terms of the power-law schemeaccounting for a HIs-network expansion and contraction within the pre- and post-inflection phase intervals and with rates of and , respectively (Fig. 5c, and see also Additional file 1: S10). On the left of the interval , both, the range and intensity of HIIS axial field component maximize as the HIs-network configuration exceeds its critical size marking the transition from the pre-inflection phase toward the post-inflection phase [69]. The energy levels associated with this phase transition are given bywhere can be replaced with its best-fitted Richards model, (see caption of Fig. 5 for Richards model parameters), providing with an estimation of the atom-packing energy (AE) (Fig. 5c). Similarly to the NaVAb case [69], AE maximization occurs in the vicinity of the narrow interval so that energetic coupling of the PMs with the VSs is dictated by the phase transition (Fig. 5c, d).

Equation r1 indicates that interatomic HIs law is robust to microscopic modifications of the atomic structure, e.g., addition, removal or deletion of a small number of atoms occurring due to a mutation-induced perturbation of .1 This happens however at the cost of re-tuning HIcIS and, hence, also AE, that in the case of small-amplitude perturbations of , are expected to be up- and down-regulated toward and away from the critical point, respectively, in a power-law fashion described by r2 (Fig. 5c). Mutation-induced perturbations propagating throughout the structure are thus expected to be amplified in the vicinity of the critical point while, on the other hand, to be damped out toward the interior (i.e., toward the HP and the SF) and toward channel exterior bounded by outer pore surface radius. Given that mutations occurring in the structural proximity of the SF are highly likely to have a deleterious LOF effect [94], observed damping-out mechanism might act as a shield protecting SF’s biological machinery from mutations occurring within the pre-inflection phase. On the other hand, mutations occurring in the post-inflection phase are unlikely to perturb the SF as they have to overcome a large energy barrier in order to reach channel interior. We thus hypothesize that critical clustering of pain-related mutation sites might actually reflect a trade-off between the two extremes; a destructive destabilization and an insignificant one.

We tested the critical-clustering hypothesis by calculating the distance of each mutation site from SF’s critical point (Additional file 1: S9b) and feeding retrieved distances into a binary classifier. We achieved to classify correctly 28 (out of 36) and 39 (out of 48) of pain-related and neutral mutation sites correctly with a cut-off distance of 5.8 Å . This translates to an area under receiver operating characteristics (ROC) curve of 0.824 and pain phenotype prediction with specificity of 0.812 and sensitivity of 0.777 (Fig. 6a). Intuitive geometrical depiction of this result requires to think of a “hot” spherical shell squeezed between the spheres of radii Å  and Å  centered at incorporating areas “a1” and “a2” thus containing the majority of pain-related mutation sites. This tendency can be deduced from Fig. 5 where we can see that correctly-classified pain-related and neutral mutation sites tend to minimize and maximize, respectively, their distance from the critical radius . The opposite holds for misclassified mutation sites. Note however that due to the pore points offset, distances of sites from line on Fig. 5 are not equal with the distances of their structural locations from the surface of the sphere of (discrepancies are of order 3.13±4.63 Å). Misclassified pain-related mutations are I136V, W1538R, I1461T, R185H, I228M, I720K, I739V and T1596I indicating that sensitivity output is qualitatively similar to the HP-based classification attempt. On the other hand, quality of specificity differs significantly among classification attempts as critical-point distance criterion misclassified neutrals M145L, M146S, R1207K, T1210N, V1613I, D890N, K1412I, K1415I and S1419N are clustering within the “hot” spherical shell in proximity to HP’s boundary.

Fig. 6

Binary classification of missense SCN9A-gene mutation sites based on their distance from SF’s critical point. a ROC curve constructed from data of distances between mutation sites and SF’s critical-point (for construction of data set see Additional file 1: S9b). Optimal threshold value 5.8 Å corresponds to specificity and sensitivity values of 0.812 and 0.777, respectively. Area under ROC curve is 0.824. b Visualization of ROC curve data. Optimal threshold value 5.8 Å is marked with black dashed line. ROC curve is constructed in R [73] by using the pROC package [93]. Two sets of missense SCN9A-gene mutation sites are employed; a pain-related set containing IEM, PPD and SFN mutation sites, and a neutral set containing mutation sites which are not expected to associate with pain disease phenotypes (Additional file 1: S8)

Finally, in order to harvest the classification power of both predictors, we linearly combined distance metrics by calculating a weighted distance average (Additional file 1: S11). The weighted distance average achieved to classify correctly 29 (out of 36) pain-related mutations and 45 (out of 48) neutrals, i.e., sensitivity = 0.805, specificity = 0.937, area under ROC curve = 0.872 (Fig. 7a). The threshold weighted distance value is 9.6 Å and it indicates which mutation sites are found in proximity to SF’s critical point and HP’s boundary.

Fig. 7

Binary classification of missense SCN9A-gene mutation sites based on a weighted distance average. a ROC curve constructed from data of distances between mutation sites and the weighted combination of SF’s critical-point and HP’s boundary (for construction of data set see Additional file 1: S11). Optimal threshold value 9.6 Å corresponds to specificity and sensitivity values of 0.805 and 0.937, respectively. Area under ROC curve is 0.872. b Visualization of ROC curve data. Optimal threshold value 9.6 Å is marked with black dashed line. ROC curve is constructed in R [73] by using the pROC package [93]. Two sets of missense SCN9A-gene mutation sites are employed; a pain-related set containing IEM, PPD and SFN mutation sites, and a neutral set containing mutation sites which are not expected to associate with pain disease phenotypes (Additional file 1: S8)

The relatively-low sensitivity of the weighted distance average is not surprising if we consider that both classifications attempts failed in correctly classifying pain-related mutation sites found far away from the HP and from the SF; misclassified pain-related mutation sites are I136V, W1538R, I1461T, R185H, I720K, I739V and T1596I and all of them are found within the post-inflection phase with the exception of I1461T which is located within the lag domain but still far away from the HP and from the SF (Figs. 3, 5a and 7). On the other hand, misclassified neutrals R1207K, V1613I and S1419N occupy “hot” spots located in proximity to SF’s critical point and HP’s boundary (Figs. 3, 5a and 7).

Discussion and concluding remarks

Criticality hypothesis in biology aims at explaining how emergence of power-laws increases biological system’s robustness and efficiency hand-in-hand with evolution. Empirical evidence for complex biological systems operating near critical points include cases of gene expression [95], DNA sequences [96], protein structures [63-66], cell growth [97] and neuronal dynamics underlying brain activity [98]. In practice, criticality implies that system dynamics are delicately balanced between an ordered state where perturbations are damped-out and a disordered state where perturbations are amplified. Consequences of critical dynamics are associated with optimal information processing [99], enhanced network stability [100] and maximal sensitivity to external stimuli [101].

In this work, instead of trying to predict the effect of missense SCN9A-gene mutations via comparing mutant NaV1.7 structures in silico, we extracted hydropathic features of the wild-type atomic environment encoding NaV1.7’s response to mutation-induced variations. Stated differently, we hypothesized that some regions of the atomic environment around NaV1.7’s pore exhibit higher sensitivity to mutation-induced perturbations due to the long-range nature of HIs guaranteeing their stability; a hallmark of SOC is that avalanche-like perturbing effects are amplified and fast-spreading throughout critical network locations [102]. To test this hypothesis we mapped mutation structural locations on their corresponding mutation sites and probed topological and scaling hydropathic characteristics of the atomic bulk around the pore. Importantly, this is possible due to the relatively-large number of pain-related mutations providing with the opportunity of structure-based mutation statistics and, consequently, identification of densely-populated (by mutation sites) structural domains.

We employed a closed-state structural model of the NaV1.7 that is constructed in silico via homology modeling procedures based on the pre-open NaVAb template (see “Methods” section). By doing so, we sought to initiate our investigations from a well-studied and precisely-engineered protonated NaV1.7 structure that has previously provided with clinically-relevant observations. The fact that the structural model in-use represents a protonated state of the NaV1.7 molecule is crucial here as cumulative hydropathicity-property moments calculations take explicitly into account hydrogen atoms. Biophysical significance of the structural model in-use was first experimentally validated in [7] where structure-based analysis successfully predicted that the IEM-related V400M and S241T mutations are energetically coupled and, hence, both exhibiting pharmacoresponsiveness to carbamazepine. The same structural model has later also been used in [30, 58] in order to explain mutation-induced electrophysiologal alternations in relation to atomic-level NaV1.7 structural changes. The main structural difference between the model in use here and the 6J8J NaV1.7 [82] structure is that the latter contains the complete DIII-DIV intracellular linker which is a hotspot for PEPD mutations. Thus working with the 6J8J NaV1.7 structure will allow to investigate two more mutation cases, namely, the PEPD-related mutations F1462V [103] and T1464I [15]. Given however that the present scheme fails in classifying correctly the PEPD-related I1461T mutation, it is rather unlikely that it will succeed with F1462V or T1464I. The reason for this misclassification is that the IS end of the channel, where the DIII-DIV linker helix is found, is predominantly hydrophilic and far away from the SF (see Figs. 3 and 5a).

The starting point of the presented procedures was the approximation of the atomic cumulative distribution function around NaV1.7’s pore demonstrating that packing of atoms follows a sigmoid pattern. The generality of the Richards model was found to be adequate for this modeling purpose verifying the sigmoid CDF hypothesis and, consequently, revealing a biphasic spatial organization of the atomic environment around the pore dictated by the spatial transition from the PM from the VSs. We showed that the pore is lined by a HP dominating within channel interior and that HIs stabilizing atom-packing around the SF are critically tuned with respect to the local inflection points. This NaV1.7 feature is shared with its evolutionary-ancestor, namely, with the pre-open NaVAb channel, suggesting that location and nature of HIs critical tuning might be conserved from NaVChs of bacterial homomers to NaVChs of mammalian heteromers [69].

Pain-related mutations tend to occupy structural locations in proximity to HP’s boundary while maintaining a critical HIs-distance from the SF. Geometrically, this result indicates that the majority of pain-related mutations are found within a spherical shell around the SF incorporating parts of the HP. What might be the evolutionary principle underlying this non-random mutation distribution around NaV1.7’s pore? Given that the hydrophilic DEKA SF sequence is conserved among human and non-human NaVCh templates [82], we propose that occurrence of mutations at critical hydropathic-interactions distance from the SF might reflect an evolutionary trade-off between potentially-deleterious destabilizations occurring too close to the SF, and insignificant polymorphisms occurring far away from it. According to this rationale, mutations occupying critical hydropathic-interactions network locations lead to a GOF effect by increasing channel’s configurational space and, consequently, expanding physiological range of ion currents, while not risking structure deletion or severe destabilizations that can induce a LOF effect [94].

Misclassification of seven pain-related events found within the post-inflection phase (namely, of I136V, W1538R, I1461T, R185H, I720K, I739V and T1596I) suggests that the destabilizing mutation effect within the post-inflection phase and, specifically, within the VSs needs to be locally investigated. In particular, misclassified pain-related events are likely to perturb local features of the VSs which are however crucial for physiological gating behavior. It might therefore be useful for future studies to consider a decoupling of the PM from the VSs in order to focus solely on the cumulative hydropathic topology and HIs-networking within the VSs. Alternatively, considering biophysical characteristics of substituted amino acids (e.g., size, charge, degree of conservation) might also contribute in improving classification accuracy as it would provide a more detailed picture of the mutation effect.

Admittedly, a limitation of this study is the small (from a statistics point of view) number of available mutation sites. To resolve this issue and provide with stronger statistical validation, we may consider in future studies to increase number of neutral and pain-related mutations, for example, by introducing NaV1.7 variants found in the Genome Aggregation Database (gnomAD) [104]. A methodological weakness is that we neglected radial hydropathic effects. In particular, even if the amplitude of the HIIS radial field component is decreasing (in comparison to the amplitude of the HIIS axial field component) for increasing molecular scale, its contribution cannot be neglected for interactions between penetrating ions species and pore walls.

In summary, our findings suggest that pathophysiological evaluation of mutation sites with respect to cumulative NaV1.7 hydropathic properties can be performed with negligible computational effort and similar or even higher accuracy to [59] (reported accuracy: 0.81) but also to the more recent study of [61] where a MLE computational pipeline was employed (reported accuracy on the human NaV1.7 template: 63.5%). Given that the formation of a narrow and hydrophilic SF followed by a wide and hydrophobic CC is a widely conserved pore-architectural motif, our observations might be relevant for other voltage-gated channel species as well. Crucially, hydropathicity-property has been previously recognized as a key-marker for predicting functional outcome of genetic defects not only in NaVChs, but also in voltage-gated calcium [94] and potassium channels [105]. Finally, in an era where MLE pipelines become increasingly popular, the phenomenological framework curated in this study could provide a phenomenological basis for biophysical interpretation of MLE-retrieved predictions regarding NaVChs pathophysiological characterization and help in tracing them back on the atomic structure in a site-specific manner as recently attempted in [106].

Additional file 1: Supplementary information.