Skip to main content

Advertisement

Network topology of NaV1.7 mutations in sodium channel-related painful disorders

Abstract

Background

Gain-of-function mutations in SCN9A gene that encodes the voltage-gated sodium channel NaV1.7 have been associated with a wide spectrum of painful syndromes in humans including inherited erythromelalgia, paroxysmal extreme pain disorder and small fibre neuropathy. These mutations change the biophysical properties of NaV1.7 channels leading to hyperexcitability of dorsal root ganglion nociceptors and pain symptoms. There is a need for better understanding of how gain-of-function mutations alter the atomic structure of Nav1.7.

Results

We used homology modeling to build an atomic model of NaV1.7 and a network-based theoretical approach, which can predict interatomic interactions and connectivity arrangements, to investigate how pain-related NaV1.7 mutations may alter specific interatomic bonds and cause connectivity rearrangement, compared to benign variants and polymorphisms. For each amino acid substitution, we calculated the topological parameters betweenness centrality (B ct ), degree (D), clustering coefficient (CC ct ), closeness (C ct ), and eccentricity (E ct ), and calculated their variation (Δ value  = mutant value -WT value ). Pathogenic NaV1.7 mutations showed significantly higher variation of |ΔB ct | compared to benign variants and polymorphisms. Using the cut-off value ±0.26 calculated by receiver operating curve analysis, we found that ΔB ct correctly differentiated pathogenic NaV1.7 mutations from variants not causing biophysical abnormalities (nABN) and homologous SNPs (hSNPs) with 76% sensitivity and 83% specificity.

Conclusions

Our in-silico analyses predict that pain-related pathogenic NaV1.7 mutations may affect the network topological properties of the protein and suggest |ΔB ct | value as a potential in-silico marker.

Background

SCN9A gene encodes the alpha-subunit of voltage-gated sodium channel NaV1.7 that is expressed in dorsal root ganglion (DRG) nociceptors and in sympathetic neurons. NaV1.7 is folded into four homologous domains, each containing six transmembrane helices (S1-S6). S1–S4 helices form the voltage-sensing domain (VSD) and highly conserved basic residues in S4 sense the electric field across the membrane. S5–S6 helices with the re-entrant extracellular loop in between form the pore domain (PD) [1]. Membrane depolarisation induces a conformational change in the VSD that, through the S4-S5 linker, is transmitted to the PD and prompt the gate to open, allowing the passage of sodium ions through the pore [2]. Opening and closing of the channel modulate the subthreshold membrane potential of nociceptors and play a key role in regulating their firing.

Missense mutations in SCN9A have been associated to a spectrum of painful conditions in humans [3], including inherited erythromelalgia (IEM), [413], paroxysmal extreme pain disorder (PEPD) [1417], and small fibre neuropathy (SFN) [18, 19]. Voltage-clamp recording, performed in transfected cell lines and DRG neurons in vitro, showed that IEM-related mutations enhance the activation of NaV1.7 through a hyperpolarising shift and a slower deactivation that keeps the channel open longer once it is activated [3], thus generating a larger-than-normal inward sodium current, with greater biophysical changes at higher temperature [20]. PEPD-related NaV1.7 mutations impair channel inactivation and prolong action potentials and repetitive nociceptor firing in response to provoking stimuli, such as stretching and exposure to cold temperatures [14, 16, 21]. NaV1.7 mutations identified in SFN patients display a spectrum of electrophysiological signatures, including impaired slow inactivation, depolarised slow and fast inactivation and enhanced resurgent currents [18].

Overall, all the disease-related NaV1.7 mutations are pro-excitatory for the NaV1.7 channel, thus increasing nociceptor excitability. For those NaV1.7 mutations that have been studied by structural modelling, the gain-of-function effect stems from functionally significant changes in the biomolecular structure of NaV1.7 channel [2224]. Accordingly, gain-of-function mutations found in IEM, PEPD, and SFN patients might be expected to produce functionally significant changes in the protein structure of NaV1.7, whereas single nucleotide polymorphisms (SNPs) or variants not associated with disease would not be expected to modify the NaV1.7 protein structure in functionally significant ways. Previous NaV1.7 structural modelling, combined with functional studies, showed that the disruption of the hydrophobic ring by the F1449V [24] or the in-frame deletion Leu955Del [22] contribute to destabilizing the NaV1.7 closed-state. These studies suggest that homology modelling is a useful tool to predict functional changes in the biomolecular structure of Nav1.7. However, the nature and extent of interatomic bond variations in NaV1.7 protein structure caused by amino acid changes have not been examined over a spectrum of mutations and SNPs.

Structural modelling combined with network theory has been widely exploited in studying protein structure to identify the emergent features of global connectivity. Indeed, several studies have used network theory to provide important insights in the local topology of interactions from a global prospective with examples from the field of allosteric communication pathways [25], protein-protein interactions [26], catalytic site residues in enzymes [27] and protein-folding mechanisms [28]. Several methods have been proposed in the literature to transform the protein structures into a network by considering: (a) the C-alpha/C-beta atoms in the amino acid residues, as in a protein backbone network [29] (b) description of the atomic contacts between residues that also feature correlated motions [3032] or (c) weak and strong non-covalent protein structure network considering atom-atom interaction at the side chain level which has been proven to provide valuable biological insights [30, 33, 34]. These studies have shown that network analysis of a protein can yield a useful method to characterize the topology of the constituent amino acid residues. Protein topologies and interaction connectivity could often produce distinct small-world networks proprieties [28, 35, 36], thus having high local connectivity of residue nodes with a smaller number of long-range residue-residue interactions.

In the present study, we aimed at elucidating specific interatomic bond variations caused by amino acid changes in NaV1.7 structure by using a network-based method. We tested the hypothesis that mutations associated with IEM, PEPD and SFN cause specific types of interatomic bonds variation of NaV1.7 that can be quantified by a network-based theoretical approach able to reduce the complexity of the three-dimensional protein architecture to one-dimensional graphs [28].

Methods

Protocol description

The overall method is summarized in Fig. 1. Our methodology can be encapsulated in a protocol that has two main components: homology modelling and topology analysis. The main steps of the current protocol are: (A) Homology modelling of NaV1.7 WT based on the bacterial NavAb sodium channel template. (B) Energy minimization and structure refinement of the protein structure (C) In-silico mutagenesis is performed for pathogenetic and control group (nABN/hSNPs) mutations (Table 1). (D) Construction of inter-residue network based on weak and strong noncovalent interactions (E) Network centrality calculation and (F) the difference between mutated and WT mutated (Δ value  = mutant value -WT value ).

Fig. 1
figure1

NaV1.7 computational protocol overview. A NaV1.7 WT homology modelling of based on the bacterial NavAb sodium channel template. B Energy minimization and structure refinement of the protein structure with YAMBER force field and FG-MD server. C In-silico mutagenesis for pathogenetic and control group (nABN/hSNPs) mutations. D Transforming NaV1.7 structure into residue interaction graphs. The construction of inter-residue network was based on interatomic bonds (hydrophobic, hydrogen bonds, salt-bridges, cation-π and π-π stacking interactions) using the commands “ListIntAtom” and “ListIntBo” via YASARA software. The de novo network construction for each mutant and WT models is achieved considering the predicted binary interatomic bonds. E-F. Network centrality calculation and their relative variation between mutant and WT (Δ value  = mutant value -WT value )

Table 1 ΔB ct values of NaV1.7 mutations associated to IEM, SFN and PEPD

NaV1.7 homology modelling

A homology model of the closed-state pore domain of the NaV1.7 was generated using the crystal structure of the bacterial Arcobacter bultzeri NaV channel NaVAb [37] as a template with the human sequence NM_002977.3 through the MEMOIR server [38]. Gap region (269-340, DI) between template-target alignment and interdomain loop regions (416-726, DI-DII; 967-1175, DII-DIII; 1458-1498, DIII-DIV) were excluded from in-silico mutagenesis (Fig. 1A). The NaVAb template shared 28% sequence identity for DI, 24% for DII, 28% for DIII and 28% for DIV (overall 27% sequence identity). The four homologous domains were modelled in the clockwise direction viewed from the extracellular side as previously suggested [39, 40]. Ab-initio modelling was performed to extend the S6 helices of the PD using the Iterative Threading ASSEmbly Refinement (I-TASSER) server [41]. The final model was subjected to energy minimization and model refinement using the YAMBER force field [42] and the Fragment-Guided Molecular Dynamics (FG-MD) server [43]. The NaV1.7 WT model was subjected to stereochemical analysis with RAMPAGE server (http://services.mbi.ucla.edu/). RAMPAGE provides results in a graphical form that shows the number of residues falling in favoured region, allowed region and in outlier region.

In-Silico Mutagenesis of NaV1.7 pathogenetic and control mutations

We performed in-silico mutagenesis via WT domain replacement of NaV1.7 mutations found in IEM, PEPD or SFN patients in which gain-of-function was demonstrated by cell electrophysiology assay and that do not alter the biophysical properties of the channel (nABN). To increase the number of control variants, we added missense SNPs identified between SCN9a homologous genes sharing >90% nucleotide sequence identity using the NCBI HomoloGene Database [44]. We constructed the phylogenetic tree of the multiple sequence alignment using ClustalW via neighbor joining method (Additional file 1: S1 Text; https://www.ebi.ac.uk/Tools/phylogeny/clustalw2_phylogeny/). The mutated models were further subjected to energy minimization and model refinement using the YAMBER force field [42] and the FG-MD server (Fig. 1B) [43]. Such hSNPs have previously been used in similar studies [4547]. All the mutations and SNPs are reported in Additional file 2: Table S1.

Transforming NaV1.7 structure into residue interaction graphs

NaV1.7 structures were transformed into mathematical graphs by identifying interatomic bonds between the amino acids. The amino acid residues form the nodes and inter-node contact interaction form the edges of the graph (Fig. 1D). We identified the interatomic bonds (hydrophobic, hydrogen bonds, salt-bridges, cation-π and π-π stacking interactions) between two residues i and j as long as the atom-atom distance between them was less than 5.0 Å using the commands “ListIntAtom” and “ListIntBo” via YASARA software (Yet Another Scientific Artificial Reality Application, www.yasara.org). Hydrophobic contacts between residues were considered in the following atom groups: (a) the first carbon of CH3-, -CH2- and CHC3 (b) sp2 carbons (phenolic rings). π-π stacking were considered between (a) sp2 carbons with a hydrogen and (b) carbon, nitrogen, oxygen or sulphur atoms in planar phenolic rings. Cation-π formation was considered to be a π-π contact with the difference being that one of the interaction partners is a cation. The de novo network construction for each mutant and WT models is achieved considering the predicted binary interatomic bonds identified through YASARA software.

Topological metrics and network visualization

We computed some of the most well-known network centrality measures for each mutant and WT network NaV1.7 graph using the Cytoscape plugin NetworkAnalyzer [48], namely:

Betweenness Centrality (B ct ) and edge Betweenness centrality (EB ct ): B ct [49] is defined as the fraction of shortest pathways between all pairs of nodes of the network that go through that node. Let G = (N, E) a graph, where N is the set of the nodes and E is the set of the edges. For each node n and m in N, let d (n, m) the distance between n and m. We define

$$ \mathrm{Betweenness}\ \mathrm{centrality}\ \left(\mathrm{n}\right)={\displaystyle {\sum}_{\mathrm{s}\ \ne \mathrm{n}\ \ne \mathrm{t}}}\frac{\upsigma_{\mathrm{s}\mathrm{t}\ }\left(\mathrm{n}\right)}{\upsigma_{\mathrm{s}\mathrm{t}}}, $$
(1)

where s, t N, σst (n) is the number of shortest paths from s to t that n lies on, and σst denotes the number of shortest paths from s to t. It accounts the importance of a node facilitating interactions between other nodes. For example, a node with high Bct can operate as a bridge on many shortest paths between other nodes in the network. It is a measure of how powerful a node is able to transfer (high B ct ) or interrupt (low B ct ) the spread of information on the fastest connection between two nodes. Similarly, the EB ct of an edge is the number of shortest paths between pairs of nodes that run along it. We define:

$$ \mathrm{Edge}\ \mathrm{Betweenness}\ \left(\mathrm{e}\right)={\displaystyle {\sum}_{{\mathrm{n}}_{\mathrm{i}}\ \in\ \mathrm{N}}}{\displaystyle {\sum}_{{\mathrm{n}}_{\mathrm{j}}\kern0.5em \in\ \mathrm{N}\backslash \left\{{\mathrm{n}}_{\mathrm{i}}\right\}\ }}{\displaystyle \sum }\ \frac{{\upsigma_{{\mathrm{n}}_{\mathrm{i}}}}_{{\mathrm{n}}_{\mathrm{j}}}\ \left(\mathrm{e}\right)\ }{{\upsigma_{{\mathrm{n}}_{\mathrm{i}}}}_{{\mathrm{n}}_{\mathrm{j}}}}, $$
(2)

Where N = set of nodes; E = set of edges; \( {\upsigma_{{\mathrm{n}}_{\mathrm{i}}}}_{{\mathrm{n}}_{\mathrm{j}}} \) = number of shortest paths between ni and nj; \( {\upsigma_{{\mathrm{n}}_{\mathrm{i}}}}_{{\mathrm{n}}_{\mathrm{j}}}\ \left(\mathrm{e}\right) \) = number of shortest paths between ni and nj which pass through e E;

Degree (D): D [49] of a node (k) is defined as the total number of nodes that it is directly connected to;

Clustering Coefficient (CC ct ): Clustering Coefficient [49] is a metric commonly employed to identify well-connected sub-components in network which represents the interconnectivity of neighbors of the node. It measures the degree to which nodes tend to cluster together and is defined as the fraction of triangles around a node among the total number of possible triangles. We define

$$ \mathrm{Clustering}\ \mathrm{Coefficient}\ \left(\mathrm{n}\right)=\frac{2{\mathrm{e}}_{\mathrm{n}}}{{\mathrm{k}}_{\mathrm{n}\ }\left({\mathrm{k}}_{\mathrm{n}} - 1\right)}, $$
(3)

where kn is the number of neighbors of n and en is the number of connected pairs between all neighbors of n;

Closeness centrality (C ct ): Cct is defined as the sum of the inverted distances, i.e. farness, to all other nodes in the graph. It captures the basic intuition that the closer a node is to all other nodes in terms of path length, the more important it is. Mathematically, C ct of a node n is defined as the inverse of the sum of shortest paths from n to all other nodes m in network. We define

$$ \mathrm{Closeness}\ \left(\mathrm{n}\right)=\frac{1}{\mathrm{average}\ \left(\mathrm{d}\left(\mathrm{n},\ \mathrm{m}\right)\right)} $$
(4)

Eccentricity (E ct ): Ect measures the distance between a node n and the most distance node m; if the E ct of the node n is low, this means that all other nodes are in proximity whereas a high E ct means that there is at least one node (and all its neighbors) that is far from node n. We define E ct maximum non-infinite length of a shortest path between n and another node in the network. We define

$$ \mathrm{Eccentricity}\ \left(\mathrm{n}\right)= \max \left\{\mathrm{d}\left(\mathrm{n},\ \mathrm{m}\right)\ :\ \mathrm{m}\kern0.5em \in \kern0.5em \mathrm{N}\right\} $$
(5)

Network centrality measure variation

For each network centrality measures we calculated the difference between mutant and WT values defined as Δvalue (Δvalue = mutant value – WT value). The NaV1.7 amino acid network was visualized using Cytoscape’s Organic layout, which is a force-directed layout algorithm similar to the Fruchterman-Reingold approach [50].

Statistical analysis

Statistical analyses were performed using the R statistical Package [51]. Data are indicated as mean ± SD. Statistical significance was determined by the Wilcoxon signed-ranked test (p <0.05). The receiver operating characteristics (ROC) curve was used to assess the discriminatory power of centrality measure variations between pathogenetic NaV1.7 mutations and control groups (nABN and hSNPs). The upper-angle of ROC corresponding to the best sensitivity and specificity was used to identify the best cut-off value.

Results

NaV1.7 interatomic structure graph design

We performed homology modelling to construct the tertiary structure of the closed-state NaV1.7 sodium channel (Fig. 1). We constructed the atomic model of NaV1.7 sodium channel using the MEMOIR server [38] based on the crystal structure of the bacterial Arcobacter bultzeri NaV channel NaVAb as a template with the human sequence NM_002977.3. The first four helices S1–S4 form the VSD and the last two helices S5–S6 form the PD (Fig. 2a and b). Gap region (269-340, S5-6 extracellular linker in DI) between template-target alignment and interdomain loop regions (416-726, DI-DII; 967-1175, DII-DIII; 1458-1498, DIII-DIV) were excluded from in-silico mutagenesis. The four homologous domains were modelled in the clockwise direction viewed from the extracellular side as suggested previously [39, 40]. Ab-initio modelling was performed to extend the S6 helices of the PD using the Iterative Threading ASSEmbly Refinement (I-TASSER) server [41]. The final model was subjected to energy minimization and model refinement using the YAMBER force field [42] and the Fragment-Guided Molecular Dynamics (FG-MD) server [43] (Additional file 3: NaV1.7 pdb file). The RAMPAGE results for the NaV1.7 model showed 88.5% residues in most favored region (Additional file 4: Figure S1), 9% (90 residues) in allowed region and 2.5% (25 residues) in outlier region. A good quality Ramachandran plot has over 90% residues in the most favoured regions [52] therefore Ramachandran plot of NaV.17 it is close to a good quality model (88.5% residues in most favoured regions).

Fig. 2
figure2

NaV1.7 structure and inter-atomic network features. a View of the sodium channel α-subunit from the intracellular side of the membrane NaV1.7 is folded into four repeated domains (DI–DIV); helices S1–S4 comprise the voltage-sensing domain (VSD); helices S5–S6 and their intracellular linker comprise the pore domain (PD). b Intramembrane view of the folded model of NaV1.7. c The graph shows the topology of the mutations found in patients with inherited erythromelalgia (IEM; red), paroxysmal extreme pain disorder (PEPD; green), small-fibre neuropathy (SFN; purple) and the amino acid substitution with no biophysical abnormalities (nABN) and homologous SNPs (light blue). Nodes represent the residues and edges of the interatomic bonds. Red and black edges represent high (red) or low (grey) edge betweenness centrality (EB ct ) values, respectively. Edge thickness are proportional to EB ct and reveal that a high number of shortest paths pass through few edges. *This mutation associates with clinical features of IEM and SFN. ǂThis mutation causes in vitro biophysics changes and in vivo symptoms common both to IEM and PEPD. The NaV1.7 amino acid network were visualized using Cytoscape’s Organic layout, which is a force-directed layout algorithm similar to the Fruchterman-Reingold approach

We performed in-silico mutagenesis for 18 mutations causing IEM, 6 mutations causing SFN, 6 mutations causing PEPD (Additional file 2: Table S1), 4 mutations not causing biophysical abnormalities (nABN) in the channel (N1245S: [53]; L1267V: [53]; V1428I and T920N: Waxman, Dib-Hajj and Mantegazza, unpublished observations) and 49 SNPs identified among human and homologous mammalian (hSNPs) SCN9A genes with >90% sequence identity (Additional file 5: S2 Text). All the disease-related mutations had previously been characterized by electrophysiological assays, and found to confer gain-of-function changes to the NaV1.7 channel (Additional file 2: Table S1). The WT and mutant NaV1.7 structures were transformed into undirected graphs by the identification of hydrophobic, cation-π and π-π stacking interactions and hydrogen bonds (H-bonds) among the amino acids. In the resulting graph, amino acids are the nodes and their interactions are the edges (Fig. 2c).

Analyses of the interatomic variations caused by gain-of-function mutations

Previous studies showed that gain-of-function mutations change the biophysical properties of the channel NaV1.7 [49, 1416, 18, 54] but the underlying interatomic variations are yet to be investigated. We analyzed the interatomic variations by calculating the network centrality parameters (B ct , D, CC ct , C ct , E ct ; see methods for detailed definitions) of WT and mutated residues and the value of the variation (Δ value  = mutant value - WT value , ΔB ct , ΔD, ΔCC ct , ΔC ct , ΔE ct ) associated with each gain-of-function NaV1.7 mutation, nABN and hSNP. B ct is a measure of the centrality of a node n defined as the fraction of shortest pathways between all pairs of nodes (s, t) of the network that go through that node n [49, 55]. D of a node n is defined as the total number of nodes that it is directly connected to [49, 55]. CC ct is a metric commonly employed to identify well-connected sub-components in network which represents the interconnectivity of neighbors of a node n [35, 49]. C ct is defined as the sum of the inverted distances of a node n, i.e. farness, to all other nodes in the graph. It captures the basic intuition that the closer a node is to all other nodes, the more important it is [56]. Eccentricity (E ct ) of a node n is the greatest distance from a node n to any other node m [55].

Figure 3a-e show the profile of the topological parameters B ct , D, CC ct , C ct , and E ct in WT and mutated residues. The graphs show that both gain-of-function mutations and nABN/hSNPs modify the D values (Fig. 3c and Additional file 6: Figure S2) and CC ct values (Fig. 3d and Additional file 7: Figure S3) in a wide range but without significant differences between the groups (gain-of-function mean D = 4.30 ± 5.15; nABN and hSNP mean ∆D = 2.27 ± 2.1; p > 0.05 by Wilcoxon signed-ranked test; gain-of-function mean ∆CC ct  = 0.15 ± 0.20; nABN and hSNP mean ∆CC ct  = 0.20 ± 0.25; p > 0.05 by Wilcoxon signed-ranked test). Smaller variations were observed in C ct values (Fig. 3e and Additional file 8: Figure S4) and E ct values (Fig. 3f and Additional file 9: Figure S5) without significant differences between the groups (gain-of-function mean ∆C ct  = 0.65 ± 0.94; nABN and hSNP mean ∆C ct  = 0.71 ± 1.51; p > 0.05 by Wilcoxon signed-ranked test; gain-of-function mean ∆E ct  = 1.53 ± 3.75; nABN and hSNP mean ∆E ct  = 2.05 ± 4.62; p > 0.05 by Wilcoxon signed-ranked test). Overall, ΔD, ΔCC ct , ΔC ct , ΔE ct did not differ significantly between gain-of-function mutations and nABN and hSNPs.

Fig. 3
figure3

Topological parameter profiles of NaV1.7 gain-of-function mutations and nABN and hSNPs. a The upper panel shows the B ct profile of gain-of-function mutation; the lower panel show the B ct profile of nABN and hSNPs. Squares indicates B ct values of WT amino acids, circles indicate B ct values of mutated amino acids. The graphs highlight that the difference between B ct value of mutated amino acids and B ct value of WT amino acids is higher in the cohort of gain-of-function (GF) mutations (upper panel) compared to control (Ctrl) nABN and hSNPs (lower panel). B ct values are multipled by 100. b The box plot shows the |ΔB ct | difference between gain-of-function mutations and the cohort of nABN and hSNP variants (mean gain-of-function |∆B ct | = 1.14 ± 1.40; nABN and hSNP ∆B ct  = 0.19 ± 0.28; ***p < 0.001 by Wilcoxon signed-ranked test). |∆B ct | values are multipled by 100; dark horizontal lines and the triangular symbol represent median and mean values respectively, with the box representing the 25th and 75th percentiles, the whiskers the 5th and 95th percentiles, and the dots the outliers. c The upper panel shows the Degree (D) profile of gain-of-function mutations; the lower panel show the D profile of nABN and hSNPs. Squares indicates D values of WT amino acids, circles indicate D values of mutated amino acids. The box plot shows the |ΔD| difference between gain-of-function (GF) mutations and the cohort of nABN and hSNP (Ctrl) variants (gain-of-function mean |∆D| = 4.3 ± 5.15; nABN and hSNP |∆D| = 2.37 ± 2.10; p > 0.05 by Wilcoxon signed-ranked test); dark horizontal lines and the triangular symbol represent median and mean values respectively, with the box representing the 25th and 75th percentiles, the whiskers the 5th and 95th percentiles, and the dots the outliers. d The upper panel shows the Clustering Coefficient (CC) profile of gain-of-function mutations; the lower panel show the CC profile of nABN and hSNPs. Squares indicates CC values of WT amino acids; circles indicate CC values of mutated amino acids. The box plot shows the |ΔCC ct | difference between gain-of-function (GF) mutations and the cohort of nABN and hSNP (Ctrl) variants (mean gain-of-function |∆CC ct | = 0.15 ± 0.20; nABN and hSNP |∆CC ct | = 0.20 ± 0.25; p > 0.05 by Wilcoxon signed-ranked test); dark horizontal lines and the triangular symbol represent median and mean values respectively, with the box representing the 25th and 75th percentiles, the whiskers the 5th and 95th percentiles, and the dots the outliers. e The upper panel shows the Closeness (Cct) profile of gain-of-function mutations; the lower panel show the Cct profile of nABN and hSNPs. Squares indicates Cct values of WT amino acids, circles indicate Cct values of mutated amino acids. The box plot shows the ΔC ct difference between gain-of-function (GF) mutations and the cohort of nABN and hSNP (Ctrl) variants (mean gain-of-function |∆C ct | = 0.6 ± 0.9; nABN and hSNP |∆C ct | = 0.7 ± 1.4; p > 0.05 by Wilcoxon signed-ranked test). C ct and ∆C ct values are multipled by 100; dark horizontal lines and the triangular symbol represent median and mean values respectively, with the box representing the 25th and 75th percentiles, the whiskers the 5th and 95th percentiles, and the dots the outliers. f The upper panel shows the Eccentricity (Ect) profile of gain-of-function mutations; the lower panel show the Ect profile of nABN and hSNPs. Squares indicates Ect values of WT amino acids, circles indicate Ect values of mutated amino acids. The box plot shows the |ΔE ct | difference between gain-of-function (GF) mutations and the cohort of nABN and hSNP (Ctrl) variants (mean gain-of-function |∆E ct | = 1.53 ± 3.75; nABN and hSNP |∆E ct | = 2.05 ± 4.62; p > 0.05 by Wilcoxon signed-ranked test); dark horizontal lines and the triangular symbol represent median and mean values respectively, with the box representing the 25th and 75th percentiles, the whiskers the 5th and 95th percentiles, and the dots the outliers

We next analysed B ct values and found that pathogenic NaV1.7 mutations are characterized by higher variations of ΔB ct compared with non-pathogenic mutations and polymorphisms (Fig. 3a; Table 1 and 2). Indeed, |ΔB ct | was significantly higher in gain-of-function mutations compared with nABN and hSNPs (gain-of-function mean ∆B ct  = 1.14 ± 1.40; nABN and hSNP mean ∆B ct  = 0.19 ± 0.28; p < 0.001 by Wilcoxon signed-ranked test; Fig. 3b). ΔB ct variations associated with Nav1.7 pathogenetic mutations and nABN variants are exemplified in the structural modeling shown in the Fig. 4.

Table 2 ΔB ct values of NaV1.7 nABN and hSNPs
Fig. 4
figure4

Structural modelling of NaV1.7 variants and their interatomic bonds. a The graph shows the NaV1.7 sodium channel topology and highlights the IEM associated mutation F216S and its intra-domain bond interaction (S3 and S4; depicted in red). b Upper left inset shows the intramembrane view of the NaV1.7 channel and the amino acid F216. Upper right inset shows network view of the four NaV1.7 channel domains (DI, purple; DII, green; DIII, light blue; DIV, orange); the topology of amino acid F216 is showed as grey node. Lower insets show the bonds of WT amino acid F216 (left) and mutated amino acid S216 (right). Hydrophobic bonds are showed in green solid lines. H-bonds are showed with yellow dashed lines. F216 (DI, red) interacts with V194, V195, F198, T202 (S3, DI) and L219 (S4, DI) via hydrophobic bonds. F216 interact via H-bonds with L213 (formed by F216[NH] and L213[CO]) and L219 (F216[CO] with L219[NH]) located in S4, DI. The mutation F216S (right) interrupts all the hydrophobic interactions with S3 residues and created new H-bonds (S216[NH] with A212[CO]) causing a decrease of B ct and EB ct values. c The graph shows the NaV1.7 sodium channel topology and highlights the IEM associated mutation L858H and its inter-domain bond interaction (S4-S5 and S4; depicted in red). d Upper left inset shows the intramembrane view of the NaV1.7 channel and the amino acid L858. Upper right insets show network view of the four NaV1.7 channeldomains (DI, purple; DII, green; DIII, light blue; DIV, orange); the topology of amino acid L858 is showed as grey node. Lower inset shows the bonds of WT amino acid L858 (left) and mutated amino acid H858 (right). Hydrophobic bonds are showed in green solid lines. H-bonds are indicated by yellow dashed lines. L858 residue (red, S4-S5; DII) interacts with I234 (DI; S4-S5), V861 (DII; S4-S5), N950, L951 and V947 (DII; S6) through hydrophobic bonds and through H-bonds with L862 (formed by L858[CO] and L862[NH]) located in DII; S4-S5. L858H mutation interrupts hydrophobic interaction with I234 (DI; S4-S5), V861 (DII; S4-S5), N950 and forms new H-bonds with A854 (formed by H858[NH] and A854[CO]) (DII; S4-S5) and V947 (formed by H858[NE2] and V947[CO]) (DII; S6). These changes decrease B ct value of amino acid 858 from 2.2 to 0.39. e The graph shows the NaV1.7 sodium channel topology and highlights the nABN mutation L1267V that is located in the domain DIII; S3 depicted in red. f Upper left inset show the intramembrane view of the NaV1.7 channel and the amino acid L1267. Upper right inset shows network view of the four NaV1.7 channel domains (DI, purple; DII, green; DIII, light blue; DIV, orange); the topology of amino acid L1267 is showed as grey node. Lower inset show the bonds of WT amino acid L1267 (left) and mutated amino acid V1267 (right). Hydrophobic bonds are showed in green solid lines. H-bonds are indicated by yellow dashed lines. L1267[NH] (VSD in DIII) interacts through H-bonds with V1263[CO]. V1267mutation interacts with V1263 through a hydrophobic bond. This change does not modify B ct value of the residue 1267

Figure 4a and b shows the B ct topological proprieties of the F216S mutation associated to IEM [11, 57]. In the WT protein, F216 is located in VSD (S4) of DI and is predicted to mediates hydrophobic interactions with V194, V195, F198, T202 (S3, DI) and L219 (S4, DI). F216 is also predicted to mediate two H-bonds: F216[NH] with L213[CO] and F216[CO] with L219[NH] residues (S4, DI). Upon mutation, the hydrophobic interaction between F216S (S4) and the S3 residues (DI; VSD) are interrupted. The H-bonds F216[NH] with L213[CO] are interrupted. New H-bonds between S216[NH] and A212[CO] are created. All these changes yield negative B ct variation (ΔB ct  = -1.71, Fig. 4a and b; Additional file 10: S1 YASARA; Additional file 11: S2 YASARA). L858H is another IEM-associated mutation [4, 9, 17]. In the WT protein, L858 is located in S4-S5 and is predicted to interacts with I234 (DI; S4-S5), V861 (DII; S4-S5), N950, L951 and V947 (DII; S6) through hydrophobic bonds and through H-bonds formed by L858[CO] and L862[NH]) (DII; S4-S5). L858H mutation interrupts hydrophobic interaction with I234 (DI; S4-S5), V861 (DII; S4-S5), N950 (DII; S6) and forms new H-bonds by H858[NH] and A854[CO] (DII; S4-S5) and by H858[CO] with V947[NH] (DII; S6) leading to a negative ΔB ct value (-1.85) (Fig. 4c and d; Additional file 12: S3 YASARA; Additional file 13: S4 YASARA). L1267V is an example of nABN variant that is located in the VSD of DIII which is highly conserved between human and SCN9A homologous genes (Additional file 5: S2 Text). L1267 interacts with V1263 through H-bonds formed by L1267[NH] and V1263[CO]. Upon mutation, V1267 forms new hydrophobic bond with V1263 which does not cause B ct variation (ΔB ct =0) (Fig. 4e and f; Additional file 14:S5 YASARA; Additional file 15:S6 YASARA).

Figure 5 shows the network inter-residue connectivity of the IEM-associated mutations I848T and N395K, both characterized by very high ΔB ct values. I848 is located in S4-S5 (DII) and I848T causes a significant hyperpolarising shift in activation, a slow deactivation and an increased response to small-ramp depolarisations in DRG nociceptors [4, 9, 17, 58, 59]. I848 is predicted to interact with S4-S5 (DII) and pore (DIII; S6) through I845 and F1435, which have with very high B ct values (3.4 and 6.6, respectively). Upon mutation, the interatomic bond interactions between DII (S4-S5) and DIII (pore; S6) are interrupted and therefore ΔB ct shifts to a negative value (-5.83) showing lower EBct values (Fig. 5a and b, Additional file 16: Figure S6). Conversely, N395K mutation forms interdomain hydrophobic (S4-S5; DI and DIV, Pore; DIV) and H-bonds (S4-S5; DIV and S6; DIV), leading to a positive ΔB ct (5) and higher EBct values.

Fig. 5
figure5

Network inter-residue connectivity of the IEM-associated mutations I848T and N395K. a The graph shows the NaV1.7 sodium channel topology and highlights the amino acids I848 (DI; S4-S5) and N395 (DI; S6). Inter-domain bond interaction are depicted in red for the IEM associated mutation I848T and in green for the IEM associated mutation N395K. b Upper panels show I848 and T848 networks, lower panels show N395 and K395 network. Bct and EBct evidence interatomic traffic over the network. Red-to-white color gradient of amino acids (nodes) represents Bct value (red represents high Bct and white low Bct). Red-to-black color gradient of edges (amino acid interatomic interactions) corresponds to EBct value (red represents high EBct and black low EBct). Hydrophobic bonds are showed in solid lines and H-bonds are indicated by dashed lines. I848 present high EB ct of connecting different parts of the NaV1.7 network. Upper right panels show that I848 interacts through hydrophobic interactions with S4-S5 (DII) and pore (DIII) through I845 and F1435 that are two residues having very high B ct values (3.4 and 6.6, respectively) and H-bonds with V852 and L844 (I848[CO] with V852[NH]; I848[NH] with L844[CO]). Note the difference of Bct of the upper left panel (I848Bct = 7.36) compared to the upper right panel (T848Bct = 1.52). I848T mutation interrupts the shortest paths within the network between DII (S4-S5) and DIII (pore) and therefore ΔB ct shifts to a negative value (-5.84). T848 interacts with F1435 through hydrophobic interactions and with S851 and L844 through H-bonds (T848[CO] with S851[NH]; T848[NH] with L844[CO]; T848[HG1] with L844[CO]). Lower panels show that N395 amino acid (red, S6 in pore module in DI) interacts with L1626 (S4-S5) via hydrophobic bond and via H-bonds formed by N395[CO] and A399[NH]and N395[NH]and F391[CO]. K395 mutation creates new hydrophobic bonds with V248 (S4-S5, DI), K398 (S6, DI), V1747 (S6, DIV), L1622 (S4-S5, DIV) and new H-bonds formed by K395[NZ] with N1751[CG] (S6, DIV) and K395[NZ] with A1625[CO] (S4-S5, DIV). These new bonds create a novel communication path within the network and thus increase Bct value of the residue 395 (K395Bct = 7.76 compared to the left panel N395Bct = 2.44). Edge thickness are proportional to EB ct and reveal that a high number of shortest paths pass through few edges

ΔBct distinguishes with high specificity pathogenic NaV1.7 mutations from variants not causing disease

The in-silico topological analyses described in Figures 4 and 5 was computed for all the pain disorder-related mutations (18 causing IEM, 6 causing SFN and 6 causing PEPD), and for all the 4 nABN and 49 hSNPs variants showed in Fig. 2c. The results showed that the only topological parameter that differs significantly between gain-of-function mutations and non-pathogenic amino acid changes is the |ΔB ct | value (Fig. 3). Indeed, 83% of nABN variants and hSNPs were characterized by |ΔB ct | values <0.26. The remaining 17% showed |ΔB ct | values >0.26 (42, V795I; 57, M1532V; 59, Y1537N; 63, V1565I; 68, V1613I; 73, I1399D; 67, T1590R; 81, D1662A; 82, G1674A) (Fig. 6a and Table 2). According to our NaV1.7 model structure, most of nABN and hSNPs, which are evolutionary variable, are located in VSD and P-loop domains and are predicted to be exposed to the lipid interface (Fig. 6b).

Fig. 6
figure6

Summary of ∆B ct values for all nABN and hSNP variants and gain-of-function mutations. a ∆B ct values of all the nABN and hSNP variants (light blue circles) and all the gain-of-function mutations (red circles) analysed in this study. Dashed lines indicates the cut-off value (ΔB ct  ± 0.26) that maximizes sensitivity and specificity. b Intra- and extracellular view of the NaV1.7 and locations of amino acids affected by gain-of-function mutations (red) that are linked with IEM, SFN and PEPD and control group variants (nABN and hSNPs; light blue). Localization of IEM, SFN and PEPD related mutation showed in red (red). *The M1532 residue shares the same position with the SFN-related M1532I mutation which causes in vitro biophysics changes and M1532V variant belongs to the control group.c Receiver operating curve (ROC) of gain-of-function mutations and control mutations (nABN and hSNPs) as a function of ΔB ct . Using a cut-off value of ± 0.26, ΔB ct correctly classified 44 out of 53 controls and 23 of 30 gain-of-function mutations yielding 76% sensitivity and 83% specificity. The area under the curve is 0.81 (95% Confidence Interval = 0.70 to 0.91)

Twenty-three out of 30 (77%) gain-of-function NaV1.7 mutations had |ΔB ct | > 0.26 and are located in VSD, Pore and S4-S5 of DI, DII, DIII and DIV domains (Table 1). The remaining 7 mutations (23%) had |ΔB ct | <0.26 (1, I136V; 2, R185H; 8, M1532I; 9, W1538R; 17, V1298F; 19, V1299F; 20, P1308L) (Fig. 6a and Table 1). These pathogenetic mutations with small |ΔB ct | variation are located in VSD of DI (2, R185H) and DIII (8, M1532I; 9 W1538R) or in S4-S5 linker of DIII (17, V1298F; 19, V1299F; 20, P1308L), are highly evolutionary conserved residues (Additional file 5: S2 Text) and are predicted to be exposed outside the core of the channel (exception: I136V; Fig. 6b-c).

According to these results, we hypothesized that ΔB ct might provide enough sensitivity and specificity to distinguish gain-of-function mutations from control variants. Using the cut-off value (ΔB ct  ± 0.26) that maximizes sensitivity and specificity, ΔB ct correctly classified 44 out of 53 controls variants (nABN and hSNPs) and 23 out of 30 gain-of-function mutations, yielding 76% sensitivity and 83% specificity. The area under the ROC curve analysis for the ΔB ct scores was 0.81 (Fig. 6c, 95% confidence interval CI = 0.70–0.91).

Discussion

Many phenomena can be modelled as collections of elements that interact through a complex set of connections. Network theory has become one of the most successful frameworks for studying these phenomena [60] and has led to major advances in our understanding of ecological systems [61], social and communication networks [62], brain connectivity [63] and metabolic and gene regulatory pathways in living cells [64].

Using network theory, protein structure can be described as mathematical graphs [28] that represent the interatomic connections. The topological features of amino acid residues, named nodes, can be described using centrality measures that define the reciprocal relationship in terms of connectivity and capability to influence other nodes within the network. We focused on the topological analysis of NaV1.7 gain-of-function mutations identified in patients with painful disorders. We considered a homology model of the NaV1.7 pore in the closed state and calculated the interaction of the nodes within the network through several measures of topology.

Our findings show that ΔB ct values tend to be significantly higher in NaV1.7 pain-related mutations than in control groups (nABN and hSNPs). B ct represents the influence that the shortest communication pathways have on the overall interatomic connections. Nodes with high B ct value could efficiently integrate signals (e.g. energy) and the reduction of B ct value caused by single amino acid substitutions suggests that the signalling transfer capability of the network is decreased. Conversely, the increase of B ct value suggests that a mutated node could facilitate the load transfer through the shortest communication pathways. Therefore, changes in ΔB ct reflect increased or decreased potential for connectivity of amino acid within the protein and provides numerical values about how single amino acid substitutions might act as a bottleneck for specific nodes linking different parts of the network. Previous studies of network topological parameters revealed that effective allosteric communications can be primarily provided by structurally stable residues that exhibit high B ct [65]. Therefore, B ct might provide a novel and useful tool for identifying allosteric hotspots in comparison with other centrality measures as previously suggested [25, 66].

Using the cut-off value (ΔB ct  ± 0.26) that maximizes sensitivity and specificity, our data show that ΔB ct correctly classified 44 out of 53 controls variants (nABN and hSNPs) and 23 out of 30 gain-of-function mutations, yielding 76% sensitivity and 83% specificity. The area under the ROC curve value for the ΔB ct scores was 0.81 (Fig. 6c, 95% confidence interval CI = 0.70–0.91). By contrast, our data show that none of other topological parameters (D, CC ct , C ct , and E ct ) differ significantly between controls and gain-of-function NaV1.7 mutations. Although these data suggest that the pain-related NaV1.7 gain-of-function mutations do not have significant effects on the degree of connectivity, local clustering connectivity of the neighbour nodes (i.e. their tendency to cluster together) and eccentricity (i.e. how far is each node from any other node within the network), it is important to consider that our results derive from homology modelling constructed on the closed-state pore domain of NaV1.7. A given residue may have a number of distinct interaction networks within the channel protein throughout the gating cycle, thus our modeling captures a snap shot of these interactions, and future studies are needed to further investigate interaction networks within the channel protein throughout the gating cycle.

Our NaV1.7 modeling also suggests a link between ΔB ct value and the buried or exposed nature of an amino acid substitution. Indeed, gain-of-function mutations predicted to be buried inside or close to the core of the channel have higher |ΔB ct | than the overall mean |ΔB ct | =1.14 (3, S211P; 4, F216S; 11, I234T; 5, I228M;12, S241T; 13, I848T; 15, L858H; 16, L858F; 24, N395K; 27, V872G; 30, A1746G) or the cut-off value (>0.26) (6, I739V; 25, V400M; 29, F1449V). Conversely, gain-of-function mutations predicted to face the lipid interface (exception: I136V) have lower |ΔB ct | than the overall mean |ΔB ct | (1.14) (14, G856D; 26, A863P; 28, M932L; 21, V1316A; 18, V1298D; 10, G1607R; 23, A1632E; exception: M1627K) or the cut-off value |ΔB ct | (<0.26) (17, V1298F; 19, V1299F; 20, P1308L; 2, R185H; 8, V1532I; 9, W1538R). Similarly, most of the control variants (nABN/hSNPs) predicted to face the lipid interface have low |ΔB ct | (<0.26) (Fig. 6b and c). Hence, in our NaV1.7 model, mutations predicted to be buried into the core of the channel show higher |ΔB ct | than those exposed at the interface of the membrane. This finding suggests that lipophilic interactions within the cell membrane may be disturbed by the mutations. Additional studies are required to more definitively assess the changes in lipophilic interactions that are produced by these mutations. Irrespective of the underlying mechanistic/molecular explanation, some pathogenic mutations would be missed by our method. Thus, ΔB ct should be regarded as a novel in-silico screening tool in addition to existing common predictive algorithms (e.g. Polyphen-2 [67], SIFT [68]) that could help in selecting pathogenetic mutations for functional testing.

Conclusions

Our findings demonstrate that most of the pathogenic NaV1.7 mutations identified in patients affected by severe painful disorders could be predicted, according to our homology modelling, to cause profound changes in the amino acid connectivity of the channel. Such modification may underpin the gain-of-function effects measurable in DRG nociceptors by electrophysiological assays. Based on these findings, we propose to consider Bct may therefore be a marker of pathogenic shift in the mutant channels, though prospective experimental studies will be required to validate its effectiveness and its biological meaning.

Abbreviations

B ct :

Betweenness centrality

CC ct :

Clustering coefficient

C ct :

Closeness

D:

Degree

DRG:

Dorsal root ganglion

EB ct :

Edge betweeness centrality

E ct :

Eccentricity

H-bonds:

Hydrogen bonds

hSNPs:

Homologous single nucleotide polymorphisms

IEM:

Inherited erythromelalgia

nABN:

Variants not causing biophysical abnormalities of the channel

PD:

Pore domain

PEPD:

Paroxysmal extreme pain disorder

ROC:

Receiver operating characteristics

SCN9A:

Sodium channel, voltage-gated, type IX, alpha subunit

SFN:

Painful small fibre neuropathy

VSD:

Voltage-sensing domain

References

  1. 1.

    Catterall WA. Voltage-gated sodium channels at 60: structure, function and pathophysiology. J Physiol. 2012;590:2577–89.

  2. 2.

    Chanda B, Bezanilla F. Tracking voltage-dependent conformational changes in skeletal muscle sodium channel during activation. J Gen Physiol. 2002;120:629–45.

  3. 3.

    Dib-Hajj SD, Yang Y, Black JA, Waxman SG. The Na(V)1.7 sodium channel: from molecule to man. Nat Rev Neurosci. 2013;14:49–62.

  4. 4.

    Cummins TR, Dib-Hajj SD, Waxman SG. Electrophysiological properties of mutant Nav1.7 sodium channels in a painful inherited neuropathy. J Neurosci. 2004;24:8232–6.

  5. 5.

    Han C, Rush AM, Dib-Hajj SD, Li S, Xu Z, Wang Y, et al. Sporadic onset of erythermalgia: a gain-of-function mutation in Nav1.7. Ann Neurol. 2006;59:553–8.

  6. 6.

    Harty TP, Dib-Hajj SD, Tyrrell L, Blackman R, Hisama FM, Rose JB, et al. Na(V)1.7 mutant A863P in erythromelalgia: effects of altered activation and steady-state inactivation on excitability of nociceptive dorsal root ganglion neurons. J Neurosci. 2006;26:12566–75.

  7. 7.

    Lampert A, Dib-Hajj SD, Eastman EM, Tyrrell L, Lin Z, Yang Y, et al. Erythromelalgia mutation L823R shifts activation and inactivation of threshold sodium channel Nav1.7 to hyperpolarized potentials. Biochem Biophys Res Commun. 2009;390:319–24.

  8. 8.

    Stadler T, O'Reilly AO, Lampert A. Erythromelalgia mutation Q875E Stabilizes the activated state of sodium channel Nav1.7. J Biol Chem. 2015;290(10):6316–25.

  9. 9.

    Yang Y, Wang Y, Li S, Xu Z, Li H, Ma L, et al. Mutations in SCN9A, encoding a sodium channel alpha subunit, in patients with primary erythermalgia. J Med Genet. 2004;41:171–4.

  10. 10.

    Dib-Hajj SD, Rush AM, Cummins TR, Hisama FM, Novella S, Tyrrell L, et al. Gain-of-function mutation in Nav1.7 in familial erythromelalgia induces bursting of sensory neurons. Brain. 2005;128:1847–54.

  11. 11.

    Drenth JPH, te Morsche RHM, Guillet G, Taieb A, Kirby RL, Jansen JBMJ. SCN9A mutations define primary erythermalgia as a neuropathic disorder of voltage gated sodium channels. J Invest Dermatol. 2005;124:1333–8.

  12. 12.

    Lee MJ, Yu HS, Hsieh ST, Stephenson DA, Lu CJ, Yang CC. Characterization of a familial case with primary erythromelalgia from Taiwan. J Neurol. 2007;254:210–4.

  13. 13.

    Drenth JPH, Waxman SG. Mutations in sodium-channel gene SCN9a cause a spectrum of human genetic pain disorders. J Clin Invest. 2007;117:3603–9.

  14. 14.

    Fertleman CR, Baker MD, Parker KA, Moffatt S, Elmslie FV, Abrahamsen B, et al. SCN9A Mutations in Paroxysmal Clinical Study Extreme Pain Disorder : Allelic Variants Underlie Distinct Channel Defects and Phenotypes. Neuron. 2006;52:767–74.

  15. 15.

    Jarecki BW, Sheets PL, Jackson JO, Cummins TR. Paroxysmal extreme pain disorder mutations within the D3/S4-S5 linker of Nav1.7 cause moderate destabilization of fast inactivation. J Physiol. 2008;586:4137–53.

  16. 16.

    Dib-Hajj SD, Estacion M, Jarecki BW, Tyrrell L, Fischer TZ, Lawden M, Cummins TR, Waxman SG. Paroxysmal extreme pain disorder M1627K mutation in human Nav1.7 renders DRG neurons hyperexcitable. Mol Pain. 2008;4:37.

  17. 17.

    Theile JW, Jarecki BW, Piekarz AD, Cummins TR. Nav1.7 mutations associated with paroxysmal extreme pain disorder, but not erythromelalgia, enhance Navbeta4 peptide-mediated resurgent sodium currents. J Physiol. 2011;589:597–608.

  18. 18.

    Faber CG, Hoeijmakers JGJ, Ahn HS, Cheng X, Han C, Choi JS, et al. Gain of function Naν1.7 mutations in idiopathic small fiber neuropathy. Ann Neurol. 2012;71:26–39.

  19. 19.

    Hoeijmakers JGJ, Han C, Merkies ISJ, Macala LJ, Lauria G, Gerrits MM, 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. 2012;135:345–58.

  20. 20.

    Han C, Lampert A, Rush AM, Dib-Hajj SD, Wang X, Yang Y, Waxman SG. Temperature dependence of erythromelalgia mutation L858F in sodium channel Nav1.7. Mol Pain. 2007;3:3.

  21. 21.

    Cheng X, Dib-hajj SD, Tyrrell L, Wright DA, Fischer TZ, Waxman SG. Mutations at opposite ends of the DIII / S4-S5 linker of sodium channel Na V 1.7 produce distinct pain disorders. Mol Pain. 2010;6:24.

  22. 22.

    Yang Y, Estacion M, Dib-Hajj SD, Waxman SG, System H, Haven W. Molecular architecture of a sodium channel S6 helix: radial tuning of the Nav1.7 activation gate. J Biol Chem. 2013;288:13741–7.

  23. 23.

    Yang Y, Dib-Hajj SD, Zhang J, Zhang Y, Tyrrell L, Estacion M, et al. Structural modelling and mutant cycle analysis predict pharmacoresponsiveness of a Nav1.7 mutant channel. Nat Commun. 2012;3:1186.

  24. 24.

    Lampert A, O’Reilly AO, Dib-Hajj SD, Tyrrell L, Wallace BA, Waxman SG. 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. 2008;283:24118–27.

  25. 25.

    Lee Y, Choi S, Hyeon C. Mapping the intramolecular signal transduction of G-protein coupled receptors. Proteins. 2014;82:727–43.

  26. 26.

    Gursoy A, Keskin O, Nussinov R. Topological properties of protein interaction networks from a structural perspective. Biochem Soc Trans. 2008;36:1398–403.

  27. 27.

    Liu R, Hu J. Computational Prediction of Heme-Binding Residues by Exploiting Residue Interaction Network. PLoS One. 2011;6:e25560.

  28. 28.

    Vendruscolo M, Dokholyan NV, Paci E, Karplus M. Small-world view of the amino acids that play a key role in protein folding. Phys Rev E Stat Nonlinear Soft Matter Phys. 2002;65:061910.

  29. 29.

    Greene LH, Higman VA. Uncovering network systems within protein structures. J Mol Biol. 2003;334:781–91.

  30. 30.

    Ghosh A, Vishveshwara S. A study of communication pathways in methionyl-tRNA synthetase by molecular dynamics simulations and structure network analysis. Proc Natl Acad Sci U S A. 2007;104:15711–6.

  31. 31.

    Seeber M, Cecchini M, Rao F, Settanni G, Caflisch A. Wordom: A program for efficient analysis of molecular dynamics simulations. Bioinformatics. 2007;23:2625–7.

  32. 32.

    Papaleo E, Renzetti G, Tiberti M. Mechanisms of intramolecular communication in a hyperthermophilic acylaminoacyl peptidase: a molecular dynamics investigation. PLoS One. 2012;7:e35686.

  33. 33.

    Bhattacharyya M, Ghosh A, Hansia P, Vishveshwara S. Allostery and conformational free energy changes in human tryptophanyl-tRNA synthetase from essential dynamics and structure networks. Proteins. 2010;78:506–17.

  34. 34.

    Bhattacharyya M, Vishveshwara S. Probing the allosteric mechanism in pyrrolysyl-tRNA synthetase using energy-weighted network formalism. Biochemistry. 2011;50:6225–36.

  35. 35.

    Watts DJ, Strogatz SH. Collective dynamics of “small-world” networks. Nature. 1998;393:440–2.

  36. 36.

    Daily MD, Upadhyaya TJ, Gray JJ. Contact rearrangements form coupled networks from local motions in allosteric proteins. Proteins Struct Funct Genet. 2008;71:455–66.

  37. 37.

    Payandeh J, Scheuer T, Zheng N, Catterall WA. The crystal structure of a voltage-gated sodium channel. Nature. 2011;475:353–8.

  38. 38.

    Ebejer J-P, Hill JR, Kelm S, Shi J, Deane CM. Memoir: template-based structure prediction for membrane proteins. Nucleic Acids Res. 2013;41:W379–83.

  39. 39.

    Li RA, Ennis IL, French RJ, Dudley Jr SC, Tomaselli GF, Marbán E. Clockwise domain arrangement of the sodium channel revealed by (mu)-conotoxin (GIIIA) docking orientation. J Biol Chem. 2001;276:11072–7.

  40. 40.

    Dudley Jr SC, Chang N, Hall J, Lipkind G, Fozzard HA, French RJ. mu-conotoxin GIIIA interactions with the voltage-gated Na(+) channel predict a clockwise arrangement of the domains. J Gen Physiol. 2000;116:679–90.

  41. 41.

    Zhang Y. I-TASSER server for protein 3D structure prediction. BMC Bioinformatics. 2008;9:40.

  42. 42.

    Krieger E, Joo K, Lee J, Lee J, Raman S, Thompson J, et al. Improving physical realism, stereochemistry, and side-chain accuracy in homology modeling: Four approaches that performed well in CASP8. Proteins. 2009;77 Suppl 9:114–22.

  43. 43.

    Zhang J, Liang Y, Zhang Y. Atomic-level protein structure refinement using fragment-guided molecular dynamics conformation sampling. Structure. 2011;19:1784–95.

  44. 44.

    Acland A, Agarwala R, Barrett T, Beck J, Benson DA, Bollin C, et al. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2014;44:D7–D19.

  45. 45.

    Care MA, Needham CJ, Bulpitt AJ, Westhead DR. Deleterious SNP prediction: Be mindful of your training data! Bioinformatics. 2007;23:664–72.

  46. 46.

    Sunyaev S, Ramensky V, Koch I, Lathe W, Kondrashov AS, Bork P. Prediction of deleterious human alleles. Hum Mol Genet. 2001;10:591–7.

  47. 47.

    Yue P, Moult J. Identification and analysis of deleterious human SNPs. J Mol Biol. 2006;356:1263–74.

  48. 48.

    Assenov Y, Ramírez F, Schelhorn SE, Lengauer T, Albrecht M. Computing topological parameters of biological networks. Bioinformatics. 2008;24:282–4.

  49. 49.

    Newman MEJ. Networks: An Introduction. Oxford: Oxford University Press; 2010.

  50. 50.

    Fruchterman TMJ, Reingold EM. Graph Drawing by Force-directed Placement. Software-Practice Exp. 1991;21:1129–64.

  51. 51.

    R Development Core Team. R. A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.r-project.org/. Accessed 1 June 2017.

  52. 52.

    Goh KM, Mahadi NM, Hassan O, Zaliha RN, Rahman RA, Illias RM. Molecular modeling of a predominant β-CGTase G1 and analysis of ionic interaction in CGTase. Biotechnology. 2008;7:418–29.

  53. 53.

    Brouwer BA, Merkies ISJ, Gerrits MM, Stephen G, Hoeijmakers JGJ, Faber CG. Painful neuropathies: the emerging role of sodium channelopathies. J Peripher Nerv Syst. 2014;19:53–65.

  54. 54.

    Lampert A, Eberhardt M, Waxman SG. Altered sodium channel gating as molecular basis for pain: contribution of activation, inactivation, and resurgent currents. Handb Exp Pharmacol. 2014;221:91–110.

  55. 55.

    Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: A software Environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13:2498–504.

  56. 56.

    Cohen E, Delling D, Pajor T, Werneck RF. Computing classic closeness centrality, at scale. Proc Second Ed ACM Conf Online Soc networks - COSN’14. 2014:37–50.

  57. 57.

    Choi JS, Dib-Hajj SD, Waxman SG. Inherited erythermalgia: Limb pain from an S4 charge-neutral Na channelopathy. Neurology. 2006;67:1563–7.

  58. 58.

    Wu MT, Huang PY, Yen CT, Chen CC, Lee MJ. A novel SCN9A mutation responsible for primary erythromelalgia and is resistant to the treatment of sodium channel blockers. PLoS One. 2013;8:e55212.

  59. 59.

    Han C, Dib-Hajj SD, Lin Z, Li Y, Eastman EM, Tyrrell L, et al. Early- and late-onset inherited erythromelalgia: genotype-phenotype correlation. Brain. 2009;132:1711–22.

  60. 60.

    Newman MEJ. The structure and function of complex networks. SIAM Rev. 2003;45:167–256.

  61. 61.

    Allesina S, Alonso D, Pascual M. A general model for food web structure. Science. 2008;320(5876):658–61.

  62. 62.

    Lazer D, Pentland A, Adamic L, Aral S, Barabási A-L, Brewer D, et al. Computational social science. Science. 2009;323:721–3.

  63. 63.

    Stam CJ. Modern network science of neurological disorders. Nat Rev Neurosci. 2014;15:683–95.

  64. 64.

    Jeong H, Albert R. The large-scale organization of metabolic networks. Nature. 2000;407:651–4.

  65. 65.

    Blacklock K, Verkhivker GM. Computational Modeling of Allosteric Regulation in the Hsp90 Chaperones: A Statistical Ensemble Analysis of Protein Structure Networks and Allosteric Communications. PLoS Comput Biol. 2014;10:e1003679.

  66. 66.

    Sengupta D, Kundu S. Role of long- and short-range hydrophobic, hydrophilic and charged residues contact network in protein’s structural organization. BMC Bioinformatics. 2012;13:142.

  67. 67.

    Adzhubei I, Jordan DM, Sunyaev SR. Predicting Functional Effect of Human Missense Mutations Using PolyPhen-2. Curr Protoc in Human Genetics. 2013;Volume Chapter 7:Unit7.20.

  68. 68.

    Flanagan SE, Patch AM, Ellard S. Using SIFT and PolyPhen to predict loss-of-function and gain-of-function mutations. Genet Test Mol Biomarkers. 2010;14:533–7.

  69. 69.

    Cheng X, Dib-Hajj SD, Tyrrell L, Waxman SG. Mutation I136V alters electrophysiological properties of the Na(v)1.7 channel in a family with onset of erythromelalgia in the second decade. Mol Pain. 2008;4:1.

  70. 70.

    Estacion M, Choi JS, Eastman EM, Lin Z, Li Y, Tyrrell L, Yang Y, Dib-Hajj SD, Waxman SG. Can robots patch-clamp as well as humans? Characterization of a novel sodium channel mutation. J Physiol. 2010;588:1915–27.

  71. 71.

    Takahashi K, Saitoh M, Hoshino H, Mimaki M, Yokoyama Y, Takamizawa M, et al. A case of primary erythermalgia, wintry hypothermia and encephalopathy. Neuropediatrics. 2007;38:157–9.

  72. 72.

    Cregg R, Laguda B, Werdehausen R, Cox JJ, Linley JE, Ramirez JD, et al. Novel mutations mapping to the fourth sodium channel domain of Nav1.7 result in variable clinical manifestations of primary erythromelalgia. Neuromolecular Med. 2013;15:265–78.

  73. 73.

    Ahn HS, Dib-Hajj SD, Cox JJ, Tyrrell L, Elmslie FV, Clarke AC, et al. A new Nav1.7 sodium channel mutation I234T in a child with severe pain. Eur J Pain. 2010;14:944–50.

  74. 74.

    Lampert A, Dib-Hajj SD, Tyrrell L, Waxman SG. Size matters: Erythromelalgia mutation S241T in Nav1.7 alters channel gating. J Biol Chem. 2006;281:36029–35.

  75. 75.

    Michiels JJ, te Morsche RH, Jansen JB, Drenth JP. Autosomal dominant erythermalgia associated with a novel mutation in the voltage-gated sodium channel alpha subunit Nav1.7. Arch Neurol. 2005;62:1587–90.

  76. 76.

    Cheng X, Dib-Hajj SD, Tyrrell L, Te Morsche RH, Drenth JPH, Waxman SG. Deletion mutation of sodium channel Na(V)1.7 in inherited erythromelalgia: enhanced slow inactivation modulates dorsal root ganglion neuron hyperexcitability. Brain. 2011;134:1972–86.

  77. 77.

    Estacion M, Yang Y, Dib-Hajj SD, Tyrrell L, Lin Z, Yang Y, et al. A new Nav1.7 mutation in an erythromelalgia patient. Biochem Biophys Res Commun. 2013;432:99–104.

  78. 78.

    Estacion M, Dib-Hajj SD, Benke PJ, Te Morsche RHM, Eastman EM, Macala LJ, et al. NaV1.7 gain-of-function mutations as a continuum: A1632E displays physiological changes associated with erythromelalgia and paroxysmal extreme pain disorder mutations and produces symptoms of both disorders. J Neurosci. 2008;28:11079–88.

  79. 79.

    Sheets PL, Jackson JO, Waxman SG, Dib-Hajj SD, Cummins TR. A Nav1.7 channel mutation associated with hereditary erythromelalgia contributes to neuronal hyperexcitability and displays reduced lidocaine sensitivity. J Physiol. 2007;581:1019–31.

  80. 80.

    Fischer TZ, Gilmore ES, Estacion M, Eastman E, Taylor S, Melanson M, et al. A novel Nav1.7 mutation producing carbamazepine-responsive erythromelalgia. Ann Neurol. 2009;65:733–41.

  81. 81.

    Choi JS, Zhang L, Dib-Hajj SD, Han C, Tyrrell L, Lin Z, et al. Mexiletine-responsive erythromelalgia due to a new Na(v)1.7 mutation showing use-dependent current fall-off. Exp Neurol. 2009;216:383–9.

  82. 82.

    Estacion M, Han C, Choi J-S, Hoeijmakers JG, Lauria G, Drenth JP, et al. Intra- and interfamily phenotypic diversity in pain syndromes associated with a gain-of-function variant of NaV1.7. Mol Pain. 2011;7:92.

  83. 83.

    Han C, Hoeijmakers JGJ, Ahn HS, Zhao P, Shah P, Lauria G, et al. Nav1.7-related small fiber neuropathy. Impaired slow-inactivation and DRG neuron hyperexcitability. Neurology. 2012;78:1635–43.

  84. 84.

    Choi JS, Boralevi F, Brissaud O, Sánchez-Martín J, Te Morsche RHM, Dib-Hajj SD, et al. Paroxysmal extreme pain disorder: a molecular lesion of peripheral neurons. Nat Rev Neurol. 2011;7:51–5.

  85. 85.

    Theile JW, Cummins TR. Inhibition of Navbeta4 peptide-mediated resurgent sodium currents in Nav1.7 channels by carbamazepine, riluzole, and anandamide. Mol Pharmacol. 2011;80:724–34.

Download references

Acknowledgements

We acknowledge the other members of the PROPANE (Probing the role of sodium channels in painful neuropathies) study group: Michela Taiana (Italy), Margherita Marchi (Italy), Raffaella Lombardi (Italy), Daniele Cazzato (Italy), Filippo Martinelli Boneschi (Italy), Andrea Zauli (Italy), Ferdinando Clarelli (Italy), Silvia Santoro (Italy), Ignazio Lopez (Italy), Angelo Quattrini (Italy), Janneke Hoeijmakers (The Netherlands), Maurice Sopacua (The Netherlands), Bianca de Greef (The Netherlands), Hubertus Julius Maria Smeets (The Netherlands), Rowida Al Momani (The Netherlands), Jo Michel Vanoevelen (The Netherlands), Ivo Eijkenboom (The Netherlands), Sandrine Cestèle (France), Oana Chever (France), Rayaz Malik (United Kingdom), Mitra Tavakoli (United Kingdom), Dan Ziegler (Germany).

Funding

The study was financed by institutional funding (IRCCS Foundation Carlo Besta Neurological Institute, Ricerca Corrente), the Italian Ministry of Health (Giovani Ricercatori, n. GR-2010/208) and the European Union 7th Framework Programme (grant 602273).

Availability of data and materials

The supporting data of this article are included within the article and its Additional files 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 and 15.

Authors’ contributions

DK perfomed designed the study, performed homology modeling, network analysis and wrote the manuscript. JS and MNX wrote and participated in data interpretation. BG performed statistical analysis and prepared the figures. YY, RLW, RS, PL, CGF, MG, ISJM, SDD, MM, SGW revised the manuscript. GL conceived the study and participated in data interpretation. All the authors participated, read and approved the manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Author information

Correspondence to Dimos Kapetis.

Additional files

Additional file 1: S1.

Text. Phylogenenetic tree of human SCN9A and homologous genes. (DOCX 35 kb)

Additional file 2: Table S1.

NaV1.7 mutations associated to IEM, SFN and PEPD. (DOCX 50 kb)

Additional file 3:

NaV1.7 pdf File. (PDB 1288 kb)

Additional file 4: Figure S1.

Ramachandran plot of NaV1.7 WT. (DOCX 233 kb)

Additional file 5: S2.

Text. Pairwise sequence alignment between human SCN9A and homologous genes. (DOCX 34 kb)

Additional file 6: Figure S2.

Degree variation (∆D) in NaV1.7 mutations compared to WT. (DOCX 1691 kb)

Additional file 7: Figure S3.

Clustering coefficient variation (∆CC ct ) in NaV1.7 mutations compared to WT. (DOCX 2798 kb)

Additional file 8: Figure S4.

Closeness Centrality variation (∆C ct ) in NaV1.7 mutations compared to WT.A (DOCX 2801 kb)

Additional file 9: Figure S5.

Eccentricity centrality variation (∆E ct ) in NaV1.7 mutations compared to WT. (DOCX 2709 kb)

Additional file 10: S1.

YASARA. Structural modelling of F216 NaV1.7 variant and their interatomic bonds. (SCE 709 kb)

Additional file 11: S2.

YASARA. Structural modelling of S216 NaV1.7 variant and their interatomic bonds. (SCE 706 kb)

Additional file 12: S3.

YASARA. Structural modelling of L858 NaV1.7 variant and their interatomic bonds. (SCE 707 kb)

Additional file 13: S4.

YASARA. Structural modelling of H858 NaV1.7 variant and their interatomic bonds. (SCE 775 kb)

Additional file 14: S5.

YASARA. Structural modelling of L1267 NaV1.7 variant and their interatomic bonds. (SCE 707 kb)

Additional file 15: S5.

YASARA. Structural modelling of V1267 NaV1.7 variant and their interatomic bonds. (SCE 706 kb)

Additional file 16: Figure S6.

Structural modelling variants and their interatomic bonds of I848T and N395K. (DOCX 747 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Kapetis, D., Sassone, J., Yang, Y. et al. Network topology of NaV1.7 mutations in sodium channel-related painful disorders. BMC Syst Biol 11, 28 (2017). https://doi.org/10.1186/s12918-016-0382-0

Download citation

Keywords

  • Sodium channel
  • Neuropathic pain
  • Structural modeling
  • Network analysis