Our methodology of molecular modeling with the ZMM program is described elsewhere, e.g., (Bruhova and Zhorov, 2010; Korkosh et al., 2014; Tikhonov and Zhorov, 2017). Briefly, we use the AMBER force field (Weiner et al., 1986) and calculate nonbonded interactions with the distance cutoff of 9 Å and a shifting function (Brooks et al., 1985). Electrostatic interactions were calculated with the distance- and environment-dependent dielectric function (Garden and Zhorov, 2010) without the distance cutoff for interactions involving ionized groups. The models were optimized by Monte Carlo energy minimizations (MCM) (Li and Scheraga, 1987) in the space of torsional angles and bond angles of prolines (Zhorov, 1981). MCM trajectories were terminated when 2,000 minimizations did not decrease energy of the apparent global minimum.

The hNav1.5 structure with quinidine (PDB ID: 6lqa) captured PD in likely inactivated state and VSDs in the activated states (Li et al., 2021). To model PD in the activated and closed states, we used as templates the open- and closed-pore structures of the NavAb channel (PDB codes 5vb8 and 5vb2, respectively) (Lenaeus et al., 2017). The structures of Nav1.5 and NavAb were 3D aligned by minimizing root mean square deviations of C^α atoms in the P1 helices, which are the most 3D-conserved elements of P-loop channels (Tikhonov and Zhorov, 2012), from C^α atoms of sequentially matching residues in the Kv1.2/Kv2.1 channel, the first eukaryotic P-loop channel whose crystal structure was obtained at a relatively high resolution below 2.5 Å (PDB ID: 2R9R) (Long et al., 2007). The NavAb and hNav1.5 channels have rather similar folding of the transmembrane (TM) and P1 helices, but different folding of extracellular loops, especially in the PD. Therefore, we imposed distance constraints between C^α atoms in the TM and P1 helices of hNav1.5 and matching C^α atoms of the template, whereas no distance constraints were imposed to the intra- and extracellular loops. A distance constraint is a flat-bottom parabolic function that imposes the energy penalty with the force constraint of 10 kcal·mol⁻¹ Å⁻² if a C^α atom in the model deviates from the template matching atom by more than 2 Å. The 3D alignment was performed by step-wise reducing distances between the C^α atoms in the model and template. At each step, the energy was MC-minimized with all degrees of freedom allowed to vary. More methodological details on large-scale transformations of PD can be found elsewhere (Korkosh et al., 2019).

To model VSDs in the resting states, we employed as a template the crystal structure of the prokaryotic sodium channel NavAb with deactivated VSDs (PDB ID: 6p6x) (Wisedchaisri et al., 2019). Methodology is similar to that used for modeling the channel with the open and closed PD.

We designate a residue by one-letter code with superscripted residue number and superscripted prefix indicating the channel segment where the residue is located. Experimental structures and models were visualized with the PyMol Molecular Graphics System, version 0.99rc (Schrödinger, New York, NY).

Earlier, we employed several bioinformatics tools and paralogue annotation method to reclassify 74 missense VUS/NP/CIP variants of w54 residues in the hNav1.5 channel as PDVs (Tarnovskaya et al., 2020). Following this publication, eight additional VUS/NP/CIP variants of the w54 residues were reported in ClinVar (version April, 2021). Therefore, currently the PDV set includes 82 variants.

We used a recent cryo-EM structure of the hNav1.5 channel obtained in nanodiscs (Li et al., 2021) to find intersegment contacts of the w54 residues. We instructed PyMol to display 3D neighbor residues within 5 Å from each query residue of set w54. Since the cryo-EM structure lacks hydrogen atoms, the distances were measured between heavy atoms. We have selected those 3D neighbors, which in the amino acid sequence are at least 5 positions away from the query residue and have sidechains facing the query residue sidechain (C^α atoms of glycines were considered as sidechains). We further ensured that each 3D neighbor of the query residue has at least one ClinVar-reported missense variant. We found 64 contact residues, which satisfied the above criteria and 85 missense variants of these residues. If two WT residues form a contact, their variants are also considered to contact each other even though they may be farther than 5 Å from each other.

Table 1 shows PDVs and their contact variants in the PD. For each PDV we show location (channel segment), all contact variants, associated syndromes, and their clinical significance. Table 2 shows analogous data for VSDs. Among 85 contact variants, 74 are VUS/NP/CIP variants and clinical significance of 11 variants is P/LP. Four pairs of residues from set w54 are contacting each other (columns “e” in Table 1). Contacts involving channel segments that undergo significant transition upon channel gating are marked “S” in column “f”.

To facilitate recognition of disease-related contacts, we schematically represented them in colored matrices. Each matrix represents an individual domain because inter-domain contacts of PDVs are rare. A matrix of contact variants for the inactivated PD (Figure 1) has four bottom rows, which represent PDVs, their location, and clinical significance; four left columns represent analogous characteristics of contact variants of PDVs. Variants associated with BrS1, LQT3 or non-specified arrhythmia are colored pink, green or yellow, respectively. Variants associated with other diseases are not colored. Variants, which are associated with more than one disease, are shown with two colors. Each cell in the matrix body links a PDV and its contact variant.

The matrix highlighted multiple contacts where mutation of either WT residue is associated with the same syndrome. A likely cause is that a disease-related contact contributes to stabilization (or destabilization) of a particular channel state; a weakened or strengthened contact would similarly affect the channel state stability regardless of the mutated contact partner. Structural analysis of all contacts involving PDVs is beyond goals of this study. Below we propose mechanisms of dysfunction for some PDVs and their contact variants.

In our models with the closed and open PD, which are described in a later section, conformations of P-loops are rather similar to those in the cryo-EM structure. And although P-loops do undergo conformational changes, e.g., upon slow inactivation, here we consider their folding as relatively state-independent. Figure 2 shows P-loops with some contacts, which apparently stabilize their folding and a matrix of respective contacts variants. Five hydrophobic residues (^IS6M³⁹⁰, ^IS6M³⁹⁷, ^IVS5M¹⁶⁶⁸, ^IVP1M¹⁷⁰¹, and L¹⁷⁰⁴) form a tight cluster (Figures 2A,B). In variants M³⁹⁰I/V, M³⁹⁴I, M¹⁶⁶⁸T, and M¹⁷⁰¹I, methionine substitutions with smaller residues (Figure 2C) would relax the cluster and increase flexibility of the smallest involved helix, IVP1. Substitutions ^IVP1L¹⁷⁰⁴H and ^IS6I³⁹⁷T would also increase flexibility of this helix. The C-end of helix IVP1 approaches the DEKA-ring alanine in the SF region (Figure 2B), which harbors the SF gate. The outer-pore of a prokaryotic sodium channel exhibits flexibility during slow inactivation (Chatterjee et al., 2018). In the Nav1.4 channel, the P-loop of domain IV is more flexible than in other domains (Tsushima et al., 1997a; Tsushima et al., 1997b). The above mutations would increase flexibility of helix IVP1, disturb the SF gate, and facilitate the channel transition to the slow-inactivated state. This proposition is supported by the fact that most substitutions in P-loops are associated with the loss-of-function syndromes, especially with BrS1 (Figure 2C).

Several other mutations involving P-loops would also increase their flexibility. For example, in the cryo-EM structure, ^IIP1F⁸⁹² and ^IIP1F⁸⁸⁸ form a stacking pair (Figure 2D) that can stabilize the α-helix (Butterfield et al., 2002). Mutation ^IIP1F⁸⁹²L affects the inter-segment contact with ^IIIS6V¹⁴⁵¹ (Table 1) and destroys the stacking pair with likely destabilization of helix IIP1. A sequential neighbor of F⁸⁹², ^IIP1R⁸⁹³ forms several strong contacts, including a salt bridge with ^IIP2E⁹⁰¹, which stabilizes mutual disposition of helices IIP1 and IIP2 (Figure 2D). Substitution of any partner in the salt bridge is associated with BrS1 (Figure 2C) likely due to destabilization of mutual disposition of helices IIP1 and IIP2. Another strong contact of ^IIP1R⁸⁹³ involves the SF glutamate ^IIP1-P2E⁸⁹⁸, which is proposed to adopt multiple conformations in a cycle of sodium ion permeation (Zhorov, 2021). The cryo-EM structure of Nav1.5 is obtained in complex with quinidine that favors inactivated channels (Li et al., 2021). However, it is unknown whether the SF gate is captured in the open or closed state. A strong electrostatic attraction to ^IIP1R⁸⁹³ may stabilize the E⁸⁹⁸ conformation in the slow-inactivated state.

Several other contacts also stabilize folding of P-loops. An example is a knob-into-hole contact ^IIIP2I¹⁴²⁴ : ^IIIP1A¹⁴¹¹ that stabilizes the mutual disposition of helices IIIP1 and IIIP2 with the selectivity-filter lysine ^IIIP1-P2K¹⁴¹⁹ (Figure 2E). Mutation I¹⁴²⁴F affects contact with A¹⁴¹¹, as well as contacts with ^IIIS5-P1F¹³⁶⁰ and ^IIIS5-P1V¹⁴⁰⁰ whose variants are associated with diseases (Table 1). Another example is a salt bridge ^IVS5-P1D¹⁶⁹⁰ : ^IP2R³⁸³ (Figures 2A,F), which is absent in BrS1-associated variants ^IVS5-P1D¹⁶⁹⁰N and R³⁸³K. In a functional study, variant D¹⁶⁹⁰N reduced I_Na and shifted the steady-state activation by +7 mV (Zeng et al., 2016).

Thus, majority of mutations affecting the P1 and P2 helices are associated with the loss-of-function BrS1. Apparently these mutations distort the SF gate and increase population of slow-inactivated channels. For comparison, two residues beyond P1 and P2 helices, ^IIS5-P1M⁸⁸¹ and ^IIS5G⁸⁵⁷, form a knob-into-hole contact that would not directly affect the SF region. Both ^IIS5-P1M⁸⁸¹I and ^IIS5G⁸⁵⁷D are PDVs (Table 1). However, unlike the above-described BrS1-associated variants that affect folding of P-helices, M⁸⁸¹I is associated with the gain-of-function LQT3, while G⁸⁵⁷D is associated with BrS1.

Figures 3A–C show cytoplasmic views of the PD, respectively, in the cryo-EM structure, closed state, and open state. Space-filled hydrophobic residues at the cytoplasmic ends of S6 helices highlight different dimensions of the activation gate in these states. Figure 3D shows intra-membrane view of the full-fledged channel model with the open PD and activated VSDs. PDVs and respective contact variants, which are specific for the models with the closed PD and open PD, are presented in Tables 2, 3.

Figure 4A shows matrix of contact variants in the cytoplasmic halves of TM helices of the open PD. A hydrophobic cluster involving ^IIIS5A¹³³⁰, ^IIIS5I¹³³⁴, ^IIIS6F¹⁴⁶⁵, ^IIIS6V¹⁴⁶⁸ and ^IVS6L¹⁷⁶¹ (Figure 4B) apparently stabilizes the open PD. Substitutions of A¹³³⁰ with bigger proline or threonine, or substitution of F¹⁴⁶⁵ with more flexible leucine are associated with the gain-of-function LQT3 likely because respective variants over-stabilize the open PD. In LQT3-associated variants L¹⁷⁶¹H/F, the open state of PD would be stabilized due to π-stacking of the aromatic substitutions with F¹⁴⁶⁵.

Variant ^IVS6F¹⁷⁶⁰C is associated with BrS1. In functional studies, F¹⁷⁶⁰C impaired fast inactivation, accelerated recovery from inactivation (Carboni et al., 2005; O'Leary and Chahine, 2015), and enhanced slow inactivation (Carboni et al., 2005). Mutant channel F¹⁷⁶⁰A, which is used as a model for drug binding assays, causes a ∼ 20 mV depolarizing shift of steady-state inactivation (Carboni et al., 2005). LQT3-associated mutation A¹³³⁰P increased the window current, whereas mutation A¹³³⁰T, which is associated with LQT3, BrS1 and other diseases, increased the window current, impaired inactivation and accelerated recovery from inactivation (Smits et al., 2005b).

Figure 4C shows the H-bond between ^IIIS6N¹⁴⁶³ and ^IIS6S⁹⁴¹ that is specific for the open PD. This is an example of inter-repeat H-bonds involving asparagine residues, which are exceptionally conserved in sodium and calcium channels. Earlier such H-bonds were predicted (Tikhonov et al., 2015) and demonstrated (Du et al., 2018) to stabilize the open conformations of these channels. LQT3-associated variant S⁹⁴¹N would over-stabilize the H-bond and hence the open state of PD, whereas variant S⁹⁴¹C, which is associated with a non-specific cardiovascular phenotype, may destabilize the H-bond and thus the channel open state. Variant S⁹⁴¹N was found in an infant died from the sudden infant death syndrome, the disease associated with both Nav1.5 gain-of-function and loss-of-function (Amin et al., 2010). The mutant channel S⁹⁴¹N expressed in Xenopus oocytes demonstrated the increased late current (Schwartz et al., 2000).

Contacts between cytoplasmic halves of TM helices in the closed PD are summarized in matrix shown in Figure 4D. Some open-PD contacts, e.g., ^IVS6L¹⁷⁶¹ : ^IIIS6F¹⁴⁶⁵, are also seen in the closed PD, but the residue orientations in the two models are substantially different (cf. Figures 4B, E). A knob-into-hole contact ^IVS6M¹⁷⁶⁶ : ^IVS5G¹⁶⁶¹, which is seen only in the closed-PD model (Table 3), and hydrophobic contact ^IVS6M¹⁷⁶⁶ : ^IVS5I¹⁶⁶⁰ (Figure 4E) apparently stabilize the closed PD conformation. It would be destabilized in LQT3-associated variants M¹⁷⁶⁶L and I¹⁶⁶⁰V where large hydrophobic residues are substituted with smaller ones. In functional studies, mutation I¹⁶⁶⁰V decreased the current, which was rescued by 48 h incubation at room temperature or by mexiletine; the variant was suggested to affect the channel folding, while mexiletine or room temperature stabilized the folding (Cordeiro et al., 2006). Mutation M¹⁷⁶⁶L significantly decreased the sodium channel expression, which was partially rescued by mexiletine; however, this mutation also demonstrated a 10-fold increase of the persistent late sodium current (Valdivia et al., 2002).

A tight contact between ^IS6A⁴¹³ and the methylene group of ^IVS6N¹⁷⁷⁴ is seen only in the closed-PD model (Table 3 and Figure 4F). Substitution of either contact partner with a larger residue would destabilize the closed conformation, resulting in LQT3 (Figure 4D). Mutant channel N¹⁷⁷⁴D expressed in the CHO-K1 cells demonstrated increased peak sodium current density, late current and enhanced activation (Kato et al., 2014). Similar electrophysiological changes, as well as markedly increased duration of the action potential, were found in the cardiac myocytes derived from induced pluripotent stem cells from a boy with LQT3-associated mutation N¹⁷⁷⁴D (Hirose et al., 2020). Variant N¹⁷⁷⁴S was found in a Chinese family with BrS1 (Du et al., 2005).

In the cryo-EM structure, ^IIIS6F¹⁴⁶⁵ forms hydrophobic contacts with four residues whose three variants are associated with LQT3 (Table 1). In our closed-PD model, ^IVS6I¹⁷⁶⁸ forms a tight contact with F¹⁴⁶⁵ whose substitution F¹⁴⁶⁵L is also associated with LQT3 (Table 3). Substitutions F¹⁴⁶⁵L or I¹⁷⁶⁸V would destabilize the hydrophobic contact and thus the closed conformation of PD. LQT3-associated substitution I¹⁷⁶⁸V would weaken contacts that stabilize both the inactivated and closed-state conformation of the PD. In functional studies, mutant channel I¹⁷⁶⁸V accelerated recovery from inactivation and increased the channel availability (Rivolta et al., 2002). In another study, mutant channel I¹⁷⁶⁸V accelerated recovery from inactivation, impaired slow inactivation, but did not affect steady-state activation and inactivation (Groenewegen et al., 2003).

PDVs and their contacts in the cryo-EM structure with the activated VSDs and models with resting PDVs are shown in Tables 4, 5, respectively. Many contacts switch upon VSD transitions between the activated-state conformations, which are captured in the cryo-EM structure, and resting-state models. We summarized these contacts in four matrices (Figures 5–7) where cells are labeled “C”, “M” or “B” depending on whether respective contacts are seen only in the cryo-EM structure, only in the resting-state model, or both in the cryo-EM structure and the model. It should be noted that ionized residues in VSDs switch contacts during channel gating. Impact of mutations in such contacts on relative stability of different channel states is difficult to predict. On the other hand, respective contacts are promising targets for further studies.

Table 7 shows domain distribution of PDVs, their contact variants, and clinical significance of the latter variants. In the PD, 27 WT residues with a total of 38 PDVs are engaged in contacts with 45 WT residues, which have a total of 63 missense disease variants. In the four VSDs, 17 WT residues with a total of 33 PDVs contact 39 WT residues with a total of 55 missense disease variants. While majority of PDVs and their contact variants are seen in the cryo-EM structures, models of the open/closed PD and resting VSDs revealed an additional PDV and 34 contact variants. Thus, in the considered channel states, w54 residues with 82 PDVs form contacts with 84 WT residues that have a total of 118 disease variants. Current number of PDVs exceeds the published number (Tarnovskaya et al., 2020) because new variants of w54 residues were recently reported in ClinVar. For majority of the published PDVs clinical significance is still VUS/NP/CIP, but two variants, ^IS5L²⁵⁰V and ^IVS6M¹⁷⁶⁶L, were recently reclassified as P/LP. As many as 104 contact disease variants are categorized as VUS/NP/CIP variants. Further bioinformatics analysis is required to predict potentially damaging effect of these variants, but the fact that they are coupled with earlier predicted PDVs suggests that a large portion of the 104 variants would be reclassified as a P/LP variants.

Table 8 shows domain distribution of disease variants. Most of these are associated with BrS1, followed by LQT3 and unspecified arrhythmia. A total of 32 variants are associated with other diseases, which are specified in Tables 1–5. In the contact variant matrices respective variants are not colored. It should be noted that most of “other” diseases are associated with the Nav1.5 loss-of-function. The cases when variants associated with “other” loss-of function diseases are coupled with BrS1-associated variants provide additional examples of contacts in which mutation of either partner destabilizes the open channel.

Only eleven of 82 PDVs (13%) are not coupled with other disease variants (Table 6). One PDV (^IS6V³⁹⁶L) is involved in tight intra-PD contacts, while 10 PDVs are exposed to lipids, cytoplasm or extracellular space. Three variants (^IIS1-S2L⁷³⁶P, ^IIS1-S2H⁷³⁸P, and ^IIIS2-S3G¹²⁶²S) affect flexibility of VSD linkers. The fact that 87% of PDVs are coupled with disease variants supports bioinformatics-based predictions of PDVs (Tarnovskaya et al., 2020) and suggests that many contact variants of PDVs, which are currently annotated as VUS/NP/CIP, may be reclassified as P/LP variants.

In the absence of cryo-EM structures of the hNav1.5 channel with the open and closed pore domain and resting VSDs, we modeled these states using as templates crystal structures of prokaryotic sodium channel NavAb (see Methods). Some of the proposed intersegment contacts in these states require future validation in functional studies. Ionizable residues in VSDs are involved in multiple state-dependent electrostatic interactions. Structural interpretation of respective mutations may be ambiguous.