Giulio Alberini,Sergio Alexis Paz,Beatrice Corradi,Cameron F. Abrams,Fabio Benfenati,Luca Maragliano
Molecular Dynamics Simulations of Ion Permeation in Human Voltage-Gated Sodium Channels
The recent determination of cryo-EM structures of voltage-gated sodium (Nav) channels has revealed many details of these proteins. However, knowledge of ionic permeation through the Nav pore remains limited. In this work, we performed atomistic molecular dynamics (MD) simulations to study the structural features of various neuronal Nav channels based on homology modeling of the cryo-EM structure of the human Nav1.4 channel and, in addition, on the recently resolved configuration for Nav1.2. In particular, single Na+ permeation events during standard MD runs suggest that the ion resides in the inner part of the Nav selectivity filter (SF). On-the-fly free energy parametrization (OTFP) temperature-accelerated molecular dynamics (TAMD) was also used to calculate two-dimensional free energy surfaces (FESs) related to single/double Na+ translocation through the SF of the homology-based Nav1.2 model and the cryo-EM Nav1.2 structure, with different realizations of the DEKA filter domain. These additional simulations revealed distinct mechanisms for single and double Na+ permeation through the wild-type SF, which has a charged lysine in the DEKA ring. Moreover, the configurations of the ions in the SF corresponding to the metastable states of the FESs are specific for each SF motif. Overall, the description of these mechanisms gives us new insights into ion conduction in human Nav cryo-EM-based and cryo-EM configurations that could advance understanding of these systems and how they differ from potassium and bacterial Nav channels.
Molecular Dynamics Simulations of Ion Permeation in Human Voltage-Gated Sodium Channels

The recent determination of cryo-EM structures of voltage-gated sodium (Nav) channels has revealed many details of these proteins. However, knowledge of ionic permeation through the Nav pore remains limited. In this work, we performed atomistic molecular dynamics (MD) simulations to study the structural features of various neuronal Nav channels based on homology modeling of the cryo-EM structure of the human Nav1.4 channel and, in addition, on the recently resolved configuration for Nav1.2. In particular, single Na+ permeation events during standard MD runs suggest that the ion resides in the inner part of the Nav selectivity filter (SF). On-the-fly free energy parametrization (OTFP) temperature-accelerated molecular dynamics (TAMD) was also used to calculate two-dimensional free energy surfaces (FESs) related to single/double Na+ translocation through the SF of the homology-based Nav1.2 model and the cryo-EM Nav1.2 structure, with different realizations of the DEKA filter domain. These additional simulations revealed distinct mechanisms for single and double Na+ permeation through the wild-type SF, which has a charged lysine in the DEKA ring. Moreover, the configurations of the ions in the SF corresponding to the metastable states of the FESs are specific for each SF motif. Overall, the description of these mechanisms gives us new insights into ion conduction in human Nav cryo-EM-based and cryo-EM configurations that could advance understanding of these systems and how they differ from potassium and bacterial Nav channels.

Voltage-gated sodium (Nav) channels are integral proteins that are responsible for Na+ conduction across the membrane of various excitable cells, where they participate in essential functions including the triggering of neuronal action potential. Several Nav mutations induce dysfunctions linked to a variety of pathological conditions such as epileptic seizures and motor neuron disease and pain, making them a major target for the design of efficient therapeutic strategies. Despite their physiological importance, little is known on their structural mechanisms, also because of the limited availability of human Nav structures. The Nav channels family includes both prokaryotic (bacterial) and eukaryotic members, which are evolutionarily divergent even if they share the same main functions. In particular, the eukaryotic subfamily of human Nav proteins comprises nine homologues, named Nav1.1 to Nav1.9, all characterized by highly conserved sequences. These channels show tissue-specific distribution, with Nav1.1, Nav1.2, and Nav1.6 expressed in the central nervous system (CNS). Moreover, unlike bacterial Nav channels, which are made up of four separated subunits, the human Nav channels are formed by four different domains, also called repeats (indicated as Rep. DI–DIV), connected together to form a single long polypeptide chain known as the ion-conducting α subunit. When viewed from the extracellular side, the DI–DIV domains are assembled in a clockwise fashion, each containing six transmembrane (TM) α helices, named S1 to S6. The S1–S4 segments of each repeat form four independent voltage-sensing domains (VSDs), while the last two helices, S5 and S6, are assembled to generate the pore domain (PD). Between the S5 and S6 segments there is a motif (named P-loop) that includes the Na+ selectivity filter (SF), which ensures the selective permeation of Na+ ions and differs between bacterial and eukaryotic channels. In particular, bacterial Nav channels display a ring of four E side chains (EEEE), as evidenced by high-resolution crystallographic structures. Several independent studies used molecular dynamics (MD) simulations to investigate the structural properties of bacterial Nav proteins. According to refs (34) and (46), the main results achieved could be summarized as follows: A first group of simulations revealed information regarding the binding sites in the SF and their multi-ion occupancy. Additional microsecond-time-scale MD simulations showed that the SF symmetric configuration of the crystal structures can be broken, affecting both the permeation and the selectivity mechanisms. The number of Na+ ions that permeate the SF is also a central issue. The observation of a multi-ion process is in agreement with the recent crystallographic structures of a bacterial NavMs channel introduced in ref (44). Moreover, the X-ray crystal structure of the closed conformation of NavAe1p, another prokaryotic Nav channel, is characterized by an outer ion site in the SF that suggests multiple ion-binding sites. Mammalian Nav channels replace the EEEE motif with a DIEIIKIIIAIV ring (here, following the notation in ref (46), subscripts indicate the Nav domain associated with each residue) that includes a lysine side chain, which selectively slows the transport of K+ over Na+. The SF is formed by the four side chains of the DIEIIKIIIAIV signature at the extracellular side and, at the inner side, by the couple of backbone carbonyl oxygen atoms belonging to the two preceding residues in each repeat. Indeed, since they are not involved in hydrogen bonds with other residues, these two inner rings of carbonyls are free to interact with the permeating cations. In addition, both experimental results and Monte Carlo simulation studies suggested the pivotal role of the E and K residues in the DEKA ring. Several investigations demonstrated that the specific position of the residues in the SF is essential to preserve the correct selectivity. Of particular interest, the DEEA variant reverts the SF to an ancestral Ca2+-selective configuration. In the case of the Nav1.2 channel, the alteration of the third residue in the DEKA motif (K1422E) affects the organism with both loss-of-function and gain-of-function effects. In addition to the DEKA signature, human channels (with the exception of Nav1.7) display a well-conserved outer ring consisting of four charged residues, EIEIIDIIIDIV (again, subscripts indicate the Nav domain associated with each residue), which is supposed to increase the electrostatic attraction and conductance of extracellular cations but not selectivity. Remarkably, the residues in both the EEDD and DEKA motifs are asymmetric in their position along the pore axis and characterized by a high degree of flexibility. To cope with such structural complexity, in the following we will refer to the region including the conductive EEDD signature, the DEKA motif, and the two inner rings as the conductivity/selectivity filter (C/SF). As for ionic occupancy, although early experiments suggested that Nav channels are characterized by multiple binding sites, simultaneous binding of multiple ions is not universally accepted, leading to disagreements on the mechanism of ionic conduction. Previous computational studies examined also the role of the lysine belonging to the DEKA ring in conduction and selectivity, but they were based on models using modified bacterial structures or homology models of mammalian Nav channels from bacterial templates. In two separate articles, the authors replaced the ring in the bacterial NavRh channel with the DEKA motif and performed MD simulations, describing the screening process of the K residue, which pushes the Na+ ion toward the D/E residues. On the other hand, the MD studies of refs (87−90), based on a homology model of the rat Nav1.4 pore-only channel from the bacterial NavAb structure, showed an essential displacement of the lysine in order to allow proper ionic permeation. More recently, the results based on extended MD simulations in ref (46) suggest that the role of lysine in conduction and selectivity is strongly dependent on its protonation state. When the K residue is deprotonated (i.e., uncharged), two- and three-ion occupancies are preferred. On the contrary, when the residue is protonated (positively charged), a reduced occupancy (one or two) is observed. Progress in understanding the molecular details of Nav channels was strongly accelerated by the recent publication of various cryo-EM structures of eukaryotic Na+ channels (for an updated summary, see ref (91)). However, the use of these structures in MD simulations remains limited and debated, mostly because their global atomic resolution is considered too low. In particular, the first two eukaryotic (nonmammalian) channels obtained, NavPaS (PDB ID 5X0M, ref (63)) and EeNav1.4 (PDB ID 5XSY, ref (64)), have a resolution of 3.8 and 4.0 Å, respectively, and many side chains belonging to the highly flexible C/SF are not resolved, hindering the investigation of processes such as ion permeation, selectivity, and drug binding. The most recent publication of the first human Nav channel (Nav1.4-β1 complex), with a global resolution of 3.2 Å and a local resolution up to 2.8 Å in the SF, represents a significant improvement in terms of quality of experimental Nav structures. However, also in this case, the use in MD simulations requires caution, since the resolution is highly variable in different parts of the protein. This notwithstanding, the growing population of cryo-EM Nav structures constitutes an important platform to study their structural and functional properties, and the number of MD-based studies is constantly growing. The exploration of ionic permeation, however, remains mostly focused on bacterial structures and derived models, and the transferability of results to eukaryotic species has been investigated only in a few works, so that a detailed thermodynamic characterization of the process in mammalian channels is still missing. In this work, we aim to investigate the residence of Na+ ions in structures of neuronal Na+ channels by means of MD simulations and free energy (FE) calculations. To this aim, we used both high-quality homology-based Nav models from cryo-EM templates and experimental cryo-EM Nav structures. The application of homology modeling is justified by the high sequence identity among human Nav (hNav) channels that spans from 85% to 91% in the TM region, thus ensuring a remarkable standard of the output configurations from the comparative modeling. In particular, starting from the first cryo-EM structure of a human Nav channel (Nav1.4, PDB ID 6AGF, at 3.20 Å resolution), we obtained the homology-modeled structures of the α subunits of three neural Nav channels (Nav1.1, Nav1.2, and Nav1.6). We performed all-atom unbiased MD simulations to assess the structural stability of these structures in the wild type (WT) configuration immersed in a water-membrane environment at a 150 mM physiological concentration of NaCl. Na+ permeation events from the external space to the internal cavity were observed in all the systems. Then, focusing on the Nav1.2 model, as a paradigmatic example of neurogical Navs, we used pore-only configurations (i.e., without the VSDs) to study the Na+ permeation in the C/SF via FE surface (FES) calculations using the temperature-accelerated MD/on-the-fly free energy parametrization (TAMD/OTFP) method. Special attention was devoted to the role of SF motif in the WT and DEEA structure, considering both the charged and uncharged version of the K or E1422 residue. The investigation was replicated on the more recent cryo-EM Nav1.2 structure (PDB ID 6J8E, at 3.00 Å resolution, based on the aforementioned Nav1.4 structure). The overall conformation of the experimental Nav1.2 is substantially identical to the previous Nav1.4 template with a root-mean-square deviation (RMSD) of 0.696 Å over 981 Cα atoms, according to the analysis exposed in ref (99). Based on this observation, not surprisingly, the results of the MD simulations performed with the latter structure nicely fit with those obtained with the starting homology-based model. Our simulations describe distinct mechanisms for single and double Na+ ions during permeation through the WT filter, with the charged K1422 residue. Remarkably, in the case of the WT configuration, both standard MD runs and FE calculations highlight single Na+ permeation events, in the course of which the ion stays in a global FE minimum in the inner part of the Nav SF. These ionic configurations observed for the WT with the charged K1422 residue are specific of the WT DEKA and differ from those obtained for the other filter motifs. Our results provide a first attempt to describe the permeation mechanisms of Na+ translocation through neuronal human Nav channels using the information provided by the recent cryo-EM hNav structures. Globally, these outcomes could help in describing the complex mechanism of ion permeation of these systems.

This section is organized in two subsections. In the first one, we report the information related to the structural modeling and the standard MD simulations, including all the details of the analysis performed. In the second subsection, we introduce the details of the OTFP/TAMD simulations. A brief descrition of the method is also included.

When we started this project, no experimental structure was available for mammalian neuronal Nav1.1, -1.2, and -1.6 channels. Therefore, we first relied on homology-based modeling to build accurate structural models for these proteins. The cryo-EM structure of the human Nav1.4 α protein, in complex with the auxiliary β1 subunit, introduced in ref (65) (PDB ID 6AGF), has been used as template. This structure is characterized by a detailed description of the pore domain, including the C/SF that regulates the permeation of ions. In particular, the critical DEKA motif is located in the inner part of the cavity and was reliably resolved with a local resolution of up to ∼2.8 Å. This value is comparable to that of other membrane protein structures employed in MD simulations. As an example, the inactive crystal structure of the μ-opioid receptor (MOR), resolved at 2.8 Å (PDB ID 4DKL(101)), was used to investigate computationally ion translocation through the allosteric sites. In Table 1, we introduce the sequence of the different domains modeled for each channel. The sequence identity between neuronal hNav channels and the hNav1.4 homologue spans from ∼85% to ∼91% in the TM region. Such values are generally considered appropriate to obtain reliable homology-based structures at the atomic resolution (see for example ref (94)). As an example, a similar criterion (∼84% identity in the TM pore region) was used in ref (105) to model the hNav1.5 channel from the cryo-EM structure of the electric eel Nav1.4 channel in order to study atomistic interactions of antiarrhythmic and local anesthetic drugs with the protein itself. Our homology-based structures include the four Rep. DI–DIV TM domains, plus the intracellular link between DIII and DIV (III–IV link). The long extracellular loop located in the DI part was not included in the modeling because it is unresolved in the template. Therefore, the residues adjacent to the missing region were connected, forming one single DI polypeptide chain.

Here we report the notation used to differentiate the Nav SF configurations investigated in this work: DEK+1A, WT configuration of the DEKA ring with the positively charged lysine. DEK0A, WT configuration of the DEKA ring with uncharged lysine. DEE–1A, mutated configuration of the DEKA ring with the charged glutamate in place of the lysine. DEE0A, mutated configuration of the DEKA ring with uncharged glutamate in place of the lysine.

The α subunit of human Nav channel is formed by four TM domains linked by cytoplasmic loops. Our TM structural models of human neuronal Nav channels were constructed in two steps following a procedure similar to the one introduced in refs (106) and (107). First, we modeled each of the four TM domains separately using the advanced protein prediction method included in the SWISS-MODEL server, with the exemption of the Rep. DIII and Rep. DIV. These last domains were modeled together, including also the cytosolic III–IV link. The conservation between the sequence of the template and the domains in each of the three neural Nav channels spans from ∼85% to ∼91%. Then, at a second stage, the whole subunit for each channel was built by structurally aligning the four domains to the template in a clockwise order viewed from extracellular side, using UCSF Chimera. The resultant tetrameric model was refined again by FG-MD to eliminate unwanted interdomain steric clashes and improve the global model quality. Before MD simulations, all the systems were aligned to the membrane normal axis z using the OPM-PPM service available at opm.phar.umich.edu/ppm_server. We show the homology-based Nav1.2 model in Figure 1, including the representation of the asymmetric Nav1.2 SF. We also show the superposition of the homology-based Nav1.2 model and the template (Figure 2) to confirm the conservation of the C/SF region. The Nav1.1 and Nav1.6 SFs are also shown in Figures S1 and S2, respectively.

We performed a 500 ns long MD simulation of the human Nav1.2 (hNav1.2) homology-based model to assess the quality of the modeling protocol. CHARMM hydrogen atoms were added with CHARMM-GUI, and disulfide bonds were assigned according to the information provided in the PDB structures 5XSY, 6AGF, and 6J8E. In particular, the disulfide bond between C950 and C959, which is located inside the pore cavity, was included in all the simulations of this work. The SF filter was modeled with a DEK+1A ring. Then the Nav1.2 structure was inserted in a homogeneous lipid patch of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) molecules with explicit water molecules, and the total charge was neutralized with a 150 mM NaCl solution, obtaining ∼290,000 atoms in total. For this MD simulation, as for the other standard MD runs described in this work, we used the NAMD software (version 2.12) and the CHARMM36m and CHARMM36 parameters for the protein and lipids, respectively, together with TIP3P model for water molecules and the associated ionic parameters with NBFIX corrections. Tetragonal periodic boundary conditions (PBCs) were applied to the simulation box to remove surface effects. Long-range electrostatic interactions were calculated using the particle mesh Ewald (PME) algorithm. Electrostatic and van der Waals interactions were calculated with a cutoff of 12 Å and the application of a smoothing decay starting to take effect at 10 Å. In this specific MD run, a time step of Δt = 2 fs was employed. To ensure maximum accuracy, electrostatic and van der Waals interactions were computed at each simulation step. All covalent bonds involving hydrogen atoms were constrained using the SHAKE and SETTLE algorithms. Before production, the system was relaxed following an equilibration for ∼35 ns by extending the CHARMM-GUI equilibration protocol in order to allow proper hydration of solvent-exposed regions of the Nav pore cavity. Then the system was simulated in the NPT ensemble using the Nosé–Hoover Langevin piston method to maintain the pressure at 1 atm and a Langevin thermostat to maintain the temperature at 310 K. The oscillation period of the piston was set at 50 fs and the damping time scale at 25 fs. The Langevin thermostat was set with a damping coefficient of 1 ps–1.

The three neuronal models Nav1.1, -1.2, and -1.6 were simulated in the WT configuration. In all these standard simulations, the K residue (K1422 in the Nav1.2 channel) was studied in the charged configuration (DEK+1A), and the hydrogen mass repartitioning (HMR) method was used to build the topology of the system. HMR is a computational methodology that redistributes the atomic masses. Thanks to this scaling, a larger time step can be used in the simulations. The method was originally proposed in ref (126) and then successfully applied to protein simulations in ref (127) using a hydrogen mass of 3 amu and a Δt of 4 fs. For our simulations, CHARMM hydrogen atoms were added with CHARMM-GUI, and disulfide bonds were assigned as done for the previous system, including the disulfide bond between C950 and C959 in Nav1.2 (C959–C968 in Nav1.1 and C944–C953 in Nav1.6, respectively), which is located inside the pore cavity. Then the Nav structures were inserted in a POPC membrane and solvated with explicit water molecules and a 150 mM NaCl solution. For all these MD simulations, we used the NAMD software and the parameters of the previous simulation with standard atomic masses. Tetragonal PBCs were applied to the simulation box to remove surface effects. Long-range electrostatic interactions were calculated using the PME algorithm. Electrostatic and van der Waals interactions were calculated with the standard CHARMM cutoff of 12 Å and the application of a smoothing decay starting to take effect at 10 Å following the instructions in ref (128). In these HMR simulations, a time step of 4 fs was employed coupled to the HMR scheme. Electrostatic and van der Waals interactions were computed at each simulation step. All covalent bonds involving hydrogen atoms were constrained using the SHAKE and SETTLE algorithms. Before the production, the systems were relaxed following an equilibration for ∼35 ns. Then, they were simulated in the NPT ensemble using the Nosé–Hoover Langevin piston method to maintain the pressure at 1 atm and a Langevin thermostat to keep the temperature at 310 K. The oscillation period of the piston was set at 300 fs and the damping time scale at 150 fs. The Langevin thermostat was employed with a damping coefficient of 1 ps–1.

We used also an experimental cryo-EM structure of the human Nav1.2 channel. The channel is bound to the peptidic μ-conotoxin KIIIA pore blocker (PDB ID 6J8E). Notably, the structure includes a Na+ ion in the C/SF, as illustrated in Figure 3. Moreover, the bottom tip of the KIIIA peptide includes only one residue, K7, whose amine group interacts with E945. The overall conformation of the experimental-based Nav1.2 is substantially identical to the Nav1.4 template with an RMSD of 0.696 Å over 981 Cα atoms, according to the analysis exposed in ref (99). Based on this observation, this structure nicely fits with our homology-based Nav1.2 model. The superposition of the two structures is illustrated in Figure S3. The comparison between the cryo-EM Nav1.2 structure and our preliminary homology-based Nav1.2 model is characterized by an RMSD of 0.746 Å over the Cα atoms belonging to the TM α helices. Furthermore, the pore-only configuration of the new Nav1.2 structure differs from our pore-only model by an RMSD of 0.534 Å. To further investigate the features of this new structure, we produced two MD runs of the α subunit of this structure (after removing all the non-α-subunit atoms, including the toxin and the ion). In particular, we performed two 500 ns MD simulations with a physiological 150 mM NaCl concentration using the HMR protocol illustrated in the previous paragraph. The RCD+ loop closure tool was used to link the adjacent residues in the missing region of the long extracellular DI loop in order to form a single DI polypeptide chain, as previously done in the homology-based models.

We calculated a set of distances between atoms of opposite repeats, reported in Tables S1−S4. For each standard MD run of the whole systems (VSDs + pore) listed in Table S5, we monitored the ionic permeation events through the SF in the presence of a physiological concentration of 150 mM of NaCl. The calculation includes the volume of the C/SF. Here we report the details of the selections using the notation of the Nav1.2 residues. Similar selections were made for the Nav1.1- and Nav1.6-based systems. The selected domain is characterized by the following features: a radius of 7 Å around the pore axis; the level of the upper circular section is defined by the z coordinate of the Cα atom following the E939 residue of the EEDD domain; the level of the lower circular section is defined by the z coordinate of the Cα atom of residue 940, further reduced by 4 Å.

We analyzed the positions within the C/SF of one sodium ion and the lysine side chain of the DEKA motif. To this goal, we combined all six MD runs of the WT systems in NaCl solution, for a total of 3 μs. For each dataset, we reconstructed the histogram of values of the z coordinate of a single Na+, provided that they were found inside the C/SF. The same analysis was performed for the nitrogen atom belonging to the DEKA lysine’s side chain, which is supposed to play a role in ionic permeation through mammalian Nav channels. The aggregation of the results from various systems is justified by the fact that the homology-based Nav1.2 model and the cryo-EM Nav1.2 structure are nearly identical in the C/SF and, because this is a highly conserved region in the various Nav proteins, the identity is transferred also to the other channels we built. Therefore, in terms of this analysis, we considered the various neuronal Nav channels as equivalent systems.

To further investigate the two-ion dynamics in the C/SF, we performed a set of additional standard MD simulations of the WT cryo-EM Nav1.2 structure (pore only) with two Na+ confined in the C/SF by applying a set of half-harmonic restraints. The cryo-EM pore was prepared using the membrane builder plugin of the CHARMM-GUI web server, embedded in a POPC bilayer and solvated with the TIP3P water model and a physiological NaCl concentration of 150 mM. Aiming at the hydration of the pore cavities, a preliminary equilibration stage of 20 ns was accomplished using the input files generated by CHARMM-GUI with the same setup as the previous simulations. In the equilibrated structure at t = 20 ns, two Na+ ions were dragged into the C/SF from the solvent in order to obtain three starting configurations accounting for two ions randomly positioned in the C/SF domain. Each configuration of the system was further subjected to an additional equilibration phase of ∼40 ns. During this stage, positional restraints were introduced on the protein with a gradually softened strength, so that in the last 10 ns only the α carbons of the S5 and S6 segments were unable to move, while P-loops were maintained free throughout the sampling. The two Na+ ions located in the C/SF region were confined by the application of two κ = 100 kcal/mol half-harmonic restraints at z > 25 Å and z < −2 Å. Simultaneously, the access to the C/SF was prevented to other ions by the application of a repulsive restraint which applies in a volume positioned approximately at the center of the DEK+1A motif by adding the following biases to the MD potential:where x = {x, y, z} are the coordinates of all the atoms in the system and the summations run over the cations to be excluded from the SF, i.e., from the volume delimited by {xmin, ymin, zmin} = {−10 Å, −10 Å, −20 Å}, and {xmax, ymax, zmax} = {10 Å, 10 Å, 25 Å}. In eqs 1, a1 = −10 Å, a2 = 10 Å, b1 = −10 Å, b2 = 10 Å, c1 = −20 Å, c2 = 25 Å, and k = 600 kcal/mol. This potential introduces a prohibitive barrier for the excluded molecules to enter the SF. For each of the equilibrated configurations, four replicas were simulated (REPLICA1–4), for a total of 3 × 4 = 12 MD runs. In REPLICA4, the exclusion volume is not applied, while the confinement of the two ions along the z axis is preserved. Every system was simulated for 100 ns in the NPT ensemble with the same parameters used during the 40 ns-long equilibration stage. The cumulative sampling amounts to 1.2 μs.

The free energy (FE) of an N-atom system, with potential energy U(x) and at temperature T, is defined aswhere the integral runs over the system configurations x ∈ and z ∈ represents a point in the lower-dimensional space of collective variables (CVs) obtained via the transformation θ: → , with M < 3N. In temperature-accelerated molecular dynamics (TAMD), a set of M auxiliary variables z = z1, ···, zM are introduced and tethered to the M CVs θ = θ1, ···, θM. The combined set (x, z) is subjected to the following potential:where V(x) is the MD force field, and κ > 0 is an adjustable spring-constant-like parameter. The extended system evolves according to coupled equations, as for example:where m is the mass, η(t) is a Gaussian process with mean 0 and covariance ⟨ηα(t)ηα′(t′)⟩ = δαα′δ(t – t′), γ̅ is a friction coefficient, and 1/β̅ is an artificial temperature with 1/β̅ > 1/β, where β = 1/kBT. As shown in ref (95) (see also ref (131)), by adjusting the parameter κ so that z ∼ θ(x(t)) and the friction coefficient γ̅ so that z moves slower than x, a trajectory z(t) can be obtained which moves at the artificial temperature 1/β̅ on the FES calculated at the physical 1/β. Indeed, it has been demonstrated that under such assumptions each auxiliary variable zj is driven by the forcewhere G(z) is given by eq 2. Hence, by choosing 1/β̅ > 1/β in eq 4, the z trajectory will rapidly visit the regions where the FE is relatively low, overcoming barriers that the system would take a long time to cross at the physical temperature 1/β. In addition, the on-the-fly free energy parametrization (OTFP) method is an efficient approach to reconstruct the FES sampled by a TAMD simulation. The idea is to introduce an approximation of the FE via a set of basis functions ϕm(z):and to determine the coefficients in eq 6 by minimizing the following error function:where λ is the vector of the λm coefficients and ⟨ ⟩ indicates the average computed during the TAMD simulation. This is equivalent to solving a linear system of equations Aλ = b, whereFollowing refs (96−98) and (100), chapeau functions are used as basis set. Other FE reconstruction methods based on the minimization of the same error function have been introduced more recently.

The NAMD code was used to integrate the atomic equations of motion in the OTFP simulations, while the auxiliary equations and the reconstruction were handled by a C code via the NAMD TcL forces plugin, available at github.com/cameronabrams/otfp. The reduced pore-only configuration of the former Nav1.2 model used for the TAMD/OTFP calculations includes the following residues: Rep. DI, from L234 to A427, with the exclusion of residues from N285 to T305; Rep. DII, from T866 to L983; Rep. DIII, from F1320 to F1476; Rep. DIV, from K1641 to I1775. The protein was aligned to the membrane normal z axis using OPM-PPM. The same pore configuration was obtained from the cryo-EM Nav1.2 structure (PDB ID 6J8E). In order to use absolute coordinates as CVs for both systems, they were aligned by superimposing the cryo-EM pore on the model one. As reference, the z coordinate of the cryo-EM K1422 NZ atom after alignment is 6.8 Å. The resolved Na+ ion in the SF (at z = 12.7 Å) was removed from the PDB file. Then, both systems were inserted in a homogeneous POPC membrane and solvated with explicit water molecules, and the total charge was neutralized with a 150 mM NaCl solution. The minimum starting size of the periodic box measured ∼ [135 × 135 × 129 ] Å3, which ensured a distance larger than 20 Å between adjacent images of the protein during the simulations. For all the TAMD/OTFP simulations, we employed the CHARMM parameters already used for the standard MD simulation with original masses. Tetragonal PBCs were applied to the simulation box, and long-range electrostatic interactions were calculated using the PME algorithm. Electrostatic and van der Waals interactions were calculated with a cutoff of 12 Å and the application of a smoothing decay starting to take effect at 10 Å. A time step of 2 fs combined with the SHAKE and SETTLE algorithms was used. In order to ensure maximum accuracy, electrostatic and van der Waals interactions were computed at each simulation step. Before production, the system was equilibrated for ∼25 ns. A snapshot of the equilibrated system is illustrated in Figure 4B. During all the TAMD/OTFP runs, soft harmonic restraints were maintained on the Cα atoms of each α S5–S6 helix in order to avoid any rigid-body rotational or translational displacement of the protein. More specifically, Cα atoms of the residues from L234 to M271 and from K399 to A427 in Rep. DI, residues from T866 to Y907 and from Q956 to L983 in Rep. DII, residues from E1321 to F1359 and from L1447 to F1476 in Rep. DIII, and residues from A1643 to F1682 and from P1749 to I1775 in Rep. DIV were restrained. On the other hand, no restraints were applied to the C/SF atoms. After the preparation, the systems were simulated in the NPT ensemble using the Nosé–Hoover Langevin piston method to maintain the pressure at 1 atm, and a Langevin thermostat to maintain the temperature at 310 K. The oscillation period of the piston was set at 50 fs and the damping time scale at 25 fs. The Langevin thermostat was employed with a damping coefficient of 1 ps–1. The ionization state of the protein residues corresponds to the neutral pH (DEK+1A and DEE–1A), except for the simulations that include the modified protonation of the K/E1422 residue (DEK0A and DEE0A). All the runs where performed in a zero-transmembrane-voltage regime. All the cations not involved in the definition of the CVs were excluded from the C/SF by applying a repulsive restraint. This restraint applies in a rectangle positioned approximately at the center of the DEK(E)A motif by including the following bias to the force field expression of the potential energy:where the summation runs over the excluded cations inside the SF, i.e., inside the rectangle delimited by the coordinates {xmin, ymin, zmin} = {−6 Å, −6 Å, −4 Å}, and {xmax, ymax, zmax} = {6 Å, 6 Å, 18 Å}, which was used for the majority of the simulations, as reported in Table S6. Here, c1 = −4 Å, c2 = 18 Å, and k = 600 kcal/mol. This potential introduces a prohibitive barrier for the excluded molecules to enter the SF. Since this volume is contained within the transmembrane domain of the channel, biases on x and y are not needed to avoid the lateral permeation of additional Na+ ions. The same potential term was previously used in ref (98) to perform TAMD/OTFP simulations of a representative voltage-gated potassium (Kv) channel while keeping water molecules outside the Kv SF. The single or two ions used for defining the CV were confined inside the filter by the application of two 300 kcal/mol half-harmonic potentials at z > 17 Å and z < 4 Å. Two-dimensional (2D) FESs were recovered by combining the following CVs: z1: the z coordinate of one Na+ in the C/SF. This CV will be denoted as rz(Na1). z12: the z coordinate of the center of mass (COM) of two Na+ ions in the C/SF. This CV will be denoted as COMz(Na1, Na2). zNZ: the z coordinate of the nitrogen atom (NZ in CHARMM notation) of the K1422 side chain. This CV will be denoted as rz(K- N H3+) in the DEK+1A and rz(K- N H2) in the DEK0A, respectively. zCD: the z coordinate of the carbon δ atom (CD in CHARMM notation) of the E1422 side chain. This CV will be denoted as rz(E- C OO–) in the DEE–1A and rz(E- C OOH) in the DEE0A, respectively. In the following, we will describe two sets of TAMD/OTFP simulations. The first, denoted with OTFP1, considers z1 and either zNZ or zCD depending on the composition of the 1422 residue. The second one, denoted with OTFP2, considers z12 and either zNZ or zCD. The auxiliary temperature was 6000 K, the auxiliary friction 800 (kcal/mol) ps/Å2, and the spring constant 1500 kcal/mol Å–2, following ref (98). The OTFP1 and OTFP2 setups were produced for all the SF configurations included in this work (DEK+1A, DEK0A, DEE–1A, and DEE0A) in the presence of a 150 mM NaCl concentration using both the homology-based Nav1.2 model and the cryo-EM Nav1.2 structure. A list of the OTFP simulations is reported in Table S6.

Using the string method, we computed the minimum free energy path (MFEP) over the averaged DEK+1A surface. The initial string was generated using the union of two linear interpolations, one between (5.0, 11.4) and (3, 10) (8 points) and the other between (3.0, 10.0) and (15.0, 7.8) (50 points). The string method was implemented with a MATLAB script, using finite differences to obtain forces from the FES. A tolerance of 10–7 was used as convergence threshold for the mean deviation between two consecutive string parametrizations. The FE profile along the MFEP was obtained by using the free energy values of the MFEP points in the FES.

In order to further test the convergence of our FESs, we also performed longer simulations using the cryo-EM Nav1.2 (OTFP1 DEK+1A) with the HMR setup. Three runs (REPLICA1/2/3) based on the same protocol of the standard HMR MD simulations, using an NPT ensemble maintained by the Nosé–Hoover Langevin piston method to keep the pressure at 1 atm and a Langevin thermostat to maintain the temperature at 310 K. Here, the oscillation period of the piston was set at 300 fs and the damping time scale at 150 fs. The Langevin thermostat was employed with a damping coefficient of 1 ps–1. Two additional runs were performed (REPLICA4/5) using a Langevin damping coefficient of 0.5 ps–1 (in REPLICA4 and REPLICA5) and a number of steps between updates of the nonbonded pair list of 5 (in REPLICA5) instead of 20, as suggested by ref (128) for improved accuracy.

All the MD trajectories were visualized and analyzed using UCSF Chimera (www.cgl.ucsf.edu/chimera/), PyMOL (pymol.org/2/), NAMD-COLVAR, and VMD (www.ks.uiuc.edu/Research/vmd/) with in-house Tcl scripts. The VMD QuickSurf representation was used to represent the volume occupancy of selected atoms in the metastable conformations extracted from the TAMD/OTFP simulations. This VMD plugin computes an isosurface extracted from the expansion of a Gaussian density kernel map at each particle position and uses a 3D grid to sum all together and get the final density. Volume densities were colored according to charged atoms from DEKA and EEDD residues (red for all the negatively charged residues). The density of the K/E1422 residue is represented according to the specific system (blue for DEK+1A, cyan for DEK0A, red for DEE–1A, orange for DEE0A). Na+ density is also included (blue). Snapshots of the equilibrated full and pore-only Nav1.2 model conformations are shown in Figure 4.

This section is organized in two parts. In the first one, we report the results obtained from the standard MD simulations, in terms of both assessment of the quality of the different structures and the evaluation of the observed Na+ permeation events through the SF. In the second subsection, we describe the thermodynamic features of Na+ translocation into the Nav1.2 SF by analyzing the 2D FESs extracted from the TAMD/OTFP simulations outlined in Materials and Methods. The list of all the simulations included in this work is reported in Figure 5 and in Tables S5 and S6.

We tested the stability of all the Nav structures involved in this study via unbiased MD simulations. For each system, we calculated the backbone RMSD of the transmembrane pore domain, including the S5 and S6 helices of each repeat (Figure S4). All the RMSD profiles show plateaus at values of ∼2 or 3 Å, thus preserving the overall pore conformation. Moreover, we calculated the RMSD of the P-loops of the different structures in each simulation in order to check the stability of the SF. These additional calculations (shown in Figure S5) confirm the stability of the C/SF domain. Not surprisingly, we found no significant difference between the homology-based model and the cryo-EM structure of the Nav1.2 channel. We also calculated a set of cross distances (CDISTs), each defined between two conserved atoms or groups of atoms belonging to opposite repeats and listed in Tables S1–S4. The average values of all CDISTs, obtained from their mean in each MD replica run, are reported in Table S7. Results for the homology-based models and the cryo-EM structure are similar. Here we introduce only the CDISTs belonging to the EEDD (d1 and d2) and the DEKA (d3 and d4) motifs, defined as follows and whose values are reported in Table 2: d1A: E387CA–D1426CA; d1B: E387CD–D1426CG d2A: E945CA–D1717CA; d2B: E945CD–D1717CG d3A: D384CA–K1422CA; d3B: D384CG–K1422CE d4A: E942CA–A1714CA; d4B: E942CD–A1714CB Here atoms are defined in the CHARMM notation: CA is the Cα atom of the backbone for each residue; CD, CG, and CE are the Cδ, Cγ, and Cε atoms of glutamate, aspartate and lysine, respectively, i.e., the most external carbon atoms in the side chains of these residues. The four CDISTs show an excellent conservation along the MD runs and confirm the stability of the C/SF filter and of the canonical initial position of the P-loops, as revealed also by the RMSD values shown in Figure S5. Selections of the residues for these additional RMSD calculations are listed in Table 3.

Despite the limited sampling of our unbiased MD simulations, they can still provide a partial picture of Na+ translocation through the filter. We thus analyzed ion occupancy in the C/SF region, as defined in Materials and Methods. In all the standard MD simulations, we observed various ion permeation events that produced the ionic occupancies reported in Table 4. While all Nav systems prefer single Na+ occupancy, they also show a small percentage of two-Na+ configurations, as further confirmed by the results exposed in Table 5, where we report the average Na+ occupancy during each standard MD run. We then looked at the specific positions of the Na+ ion in the C/SF in all the frames that showed a single cation permeating the filter. These snapshots form the majority of the ionic configurations revealed by our MD runs. To this aim, we combined the single-ion configurations from all six MD runs of the different systems. The aggregation of results from the three Navs is justified by the following observations: 1. The residues forming the C/SF domain are conserved across the three isoforms. 2. The starting, equilibrated structures of Nav1.1 and Nav1.6 have an RMSD of the P-loops less than or about 0.1 Å with respect to the analogous Nav1.2 conformation. 3. The RMSD values of the P-loops of all isoforms are stable during all the standard MD simulations (as shown in Figure S5), with maximum average value of 1.95 Å and maximum variance of 0.49 Å. 4. All isoforms show the same trend in terms of ion occupation in the C/SF domain, as shown in Tables 4 and 5. In Figure 6 we show the histograms of the z coordinate of (A) the single Na+ ion in the C/SF and (B) the nitrogen of K1422 side chain, which is considered to play a relevant role in ionic permeation through mammalian Navs. In addition, we show the histograms of the z coordinate of the two Na+ ions and of their center of mass (COM) in the cases of double occupancy. The interval of the z coordinate is between 3 and 18 Å, accordingly with the sampling of our TAMD/OTFP simulations (see below). In Figure S6A–C we also report the same histograms using an extended range of z coordinates, between −2 and 22 Å. Results show that the single Na+ mostly occupies the lower portion of the SF, below the DEKA ring and in proximity of the inner rings of carbonyls. Conversely, the lysine side chain resides at a higher position along the axis, where the positive charge is expected to interact with the negatively charged residues of the ring. These interactions involving the lysine might even occlude the SF, competing with the passage of other alkali ions from the extracellular environment. A set of six movies showing different events of Na+ permeation throughout the C/SF and its residency in proximity of the lower rings are provided at https://github.com/BeatriceCorradi/Nav_Channels_movies.git. Our findings of single Na+ in proximity of the inner region of the SF are consistent with what was observed in ref (93) via MD simulations of multiple Nav channels (including 1.2, 1.4, and 1.7) starting from cryo-EM structures. Notably, however, this position is different from the one occupied by a single Na+ ion in the Nav1.2 structure (PDB ID 6J8E). However, the cryo-EM structure captures the channel with the C/SF occluded by the bound μ-conotoxin KIIIA, which limits a comparison with our computational investigations. Calculations for the z coordinates of the two Na+ ions in the case of double occupancy show the main peak of the COM at the center of the SF domain at ∼11 Å with one Na+ in proximity of the EEDD domain and the other in the internal part of the C/SF region below the DEKA motif. Because of the limited sampling of the two Na+ configurations in our standard MDs of whole channels, we performed a set of additional simulations of the pore-only configuration of the more recent cryo-EM Nav1.2 structure by confining two cations in the C/SF and calculated their positions along the z axis during a cumulative time of 1.2 μs. The resulting positional distributions confirm the peak for the COM occupancy at ∼11 Å and the associated positions of the ions, as shown in Figure S6D. The corresponding trajectories of the two ions and the lysine in the 12 simulations are reported in Figures S7–S9.

We then used TAMD/OTFP simulations to reconstruct the thermodynamics of Na+ permeation through the C/SF of the Nav1.2 channel structures with different residue configurations, namely DEK+1A, DEK0A, DEE–1A, and DEE0A (see Materials and Methods). We performed two sets of TAMD/OTFP simulations, named OTFP1 and OTFP2, that differ in the definition of the CVs and describe, respectively, single- and two-Na+ translocation, the most probable spontaneous events observed in the standard MD simulations. To investigate possible distortions in the pore-only configuration induced by the removal of the four VSDs, we also performed a standard MD simulation, with no restraints applied, of the cryo-EM pore-only structure (∼190 ns, HMR setup). The RMSD calculations shown in Figure S10 demonstrate that the deviation of the structure from its starting state is comparable to that of the whole channels (compare with Figure S4), even for the P-loops only (compare with Figure S5). These findings support the use of the pore-only configuration in our additional MD simulations with two ions confined in the C/SF (see previous section) and also in the accelerated simulations. During all the OTFP/TAMD simulations, we monitored the stability of the C/SF domain by calculating the RMSD of the backbone for the residues belonging to the P-loops (according to the selection reported in Table 3). The calculations of the DEKA systems are shown in Figures S11 and S12 for the OTFP1 and the OTFP2, respectively. The calculations of the DEEA systems are illustrated in Figures S13 and S14, again for the OTFP1 and the OTFP2, respectively. They show a maximum average RMSD of 1.32 Å and maximum variance of 0.22 Å, in line with the data obtained from the standard MD simulations, demonstrating that the setup of the FE calculation (bias, effective temperature, exclusion volumes for other ions, absence of the VSDs) does not induce any distortion in the structure of the system. In the next two paragraphs we show the average FESs obtained from multiple independent OTFP1 and OTFP2 runs, respectively, while results for the individual trajectories are reported in Figures S15–S18. All TAMD/OTFP runs are summarized in Table S6. In all simulations, the ions not involved in the CVs definition were excluded from the C/SF (see Materials and Methods). In order to investigate the effect of this setting, we performed a set of independent simulations for the DEK+1A system using a larger exclusion domain. Results are reported in Figure S19 and do not show any remarkable difference with the smaller domain. Convergence was checked by comparing FESs obtained at different time lengths during the simulations. For instance, in Figure S20 we present the evolution of the FESs during a single replica simulation for the DEK+1A system. It is possible to see that the FES at 154 ns is very similar to the FES at 115 ns, and continuing the simulation did not result in major modifications. It is important to note that while this last figure corresponds to a single replica, the statistics of three simulations for the cryo-EM structure and three replicas for the homology model were combined to get the final FESs presented in the article, all based on the use of the standard atomic masses and the time step Δt = 2 fs, in order to have maximum accuracy. We also performed five HMR OTFP1 simulations of the cryo-EM Nav1.2 structure (DEK+1A SF), each for a time length of ∼300 ns. The associated FESs, shown in Figure S21, are very similar to those with standard atomic masses and Δt = 2 fs (OTFP1 runs for the DEK+1A motif in Figure 7). Moreover, these simulations preserved the starting configurations of the P-loops, as confirmed by the RMSD calculations shown in Figure S22. Comparison of the total energy in the three replicas with standard atomic masses (Δt = 2 fs) and the five replicas with HMR (Δt = 4 fs), shown in Figure S23, demonstrates that using TAMD/OTFP and HMR does not affect the stability of the MD integrator.

For each of the DEK+1A, DEK0A, DEE–1A, and DEE0A Nav1.2 systems, we ran six independent OTFP1 calculations, three starting from the homology model (Figure S15) and three starting from the cryo-EM structure (Figure S16). Then we calculated an average map by taking, for each point in CV space where the FESs are defined, the average of their FE values. The resulting averaged maps are reported in Figure 7, along with representative snapshots for the metastable and other states. Interestingly, the DEK+1A system shows the broader low-energy region, where it is possible to identify a gradually decreasing energy path for the ion through the SF and different configurations for the side chain of K1422. Conversely, all other systems display a single isolated minimum, each one associated with a unique orientation of the side chain of residue 1422. For the DEK+1A system, the positions of the single Na+ and K1422 at the energy minima are well-correlated with the high-probability locations obtained from our standard MD simulations (see Figure 6) and with what is described in ref (93). In these configurations (see the states y1 and y2 in Figure 7), the Na+ is located at the inner portion of the SF (rz(Na1) ≈ 4 Å), in the proximity of the internal carbonyl oxygen atoms belonging to the two residues underneath the DEKA ring. The K1422 ammonium group is found either at the level of the DEKA backbone (at rz(K-NH3+) ≈ 11 Å, y1) or at a lower position (rz(K-NH3+) ≈ 8 Å, y2). Using the string method, we computed the minimum free energy path (MFEP) over the DEK+1A surface, shown as a white trace over the corresponding map in Figure 8, where we also report the FE profile along the path. A video with snapshots of the representative configuration and volume occupancies along the entire MFEP is given in the Supporting Information. The MFEP allows us to reconstruct the following mechanism for single ion permeation: Na+ access to the upper part of the C/SF is possible only when K1422 is in a lower position, interacting with the upper ring of the internal carbonyls (configurations y4 and y3 in Figure 7); during ion permeation through the filter (rz(Na1) from 15 to 5 Å), the K1422 side chain switches from the downward to the upward conformation. When K1422 starts moving, the ion crosses it and reaches lower positions, making space for the side chain to move up. Once this is in the upward configuration, the ion finally adjusts to a slightly higher position. In terms of FE values, this last part of the path is essentially a descent toward the minimum. This picture is consistent with the role of the protonated K1422 suggested in refs (70) and (146), where it is proposed that in the up position the charged lysine interacts with the negatively charged side chain of the DEKA ring to block the access of ions from the extracellular space. Opposite to the DEK+1A configuration, all the other systems’s FESs show an isolated minimum, again with the ion located in the region surrounded by the carbonyl oxygen atoms below the DE*A ring (where * indicates the variable residue among the different configurations), at 4 Å < rz(Na1) < 8 Å. In all three cases, the modification of the K1422 residue results in higher energy values in the upper part of the C/SF. The system with the DEK0A ring is characterized by a single metastable state in the FES, y5 (Figure 7), defined by the Na+ ion at rz(Na1) ≈ 6 Å and the rz(K-NH2) at ∼7 Å. Interestingly, other works investigated the role of different protonation states of the DEKA lysine in ionic permeation. While extensive unbiased MD simulations of a chimeric bacterial/eukaryotic channel suggested that the uncharged lysine plays a minor role in the selectivity process, Monte Carlo simulations on multiple eukaryotic structures proposed that the uncharged lysine is responsible for escorting the Na+ ion to the inner carbonyl rings, resulting in a configuration that is very similar to the one identified in our DEK0A TAMD/OTFP simulations (state y5). Also the K1422E mutation (either DEE–1A or DEE0A) affects the thermodynamic features of Na+ translocation. In the DEE–1A minimum (y6), the Na+ ion is localized at rz(Na1) ≈ 8 Å, with the E1422 side chain confined in the upper region of the SF, at rz(E-COO–) ≈ 14 Å. In the case of DEE0A, the FE minimum (y7) is located at rz(Na1) ≈ 7 Å and rz(E-COOH) ≈ 9 Å, i.e., with the E1422 CD atom in the proximity of the inner ring’s backbone, in a configuration similar to that of the DEK0A system (y5) but with the ion and the glutamate side chain in a slightly higher position.

We then computed 2D FESs using the OTFP2 setup of CVs illustrated in Materials and Methods, i.e., by considering as CVs the vertical positions of the COM of the two ions and of the side chain of residue 1422. Indeed, multiple Na+ ion configurations were observed in a chimeric bacterial/mammalian model including the Nav1.2 filter and extensively discussed in many works presenting experimental results. As detailed in Table S6, and similarly to what was described for the OTFP1 simulations, for each (DEK+1A, DEK0A, DEE–1A, DEE0A) system we performed six independent OTFP2 calculations, three using the homology model (Figure S17) and three the cryo-EM structure (Figure S18), and averaged the results in a single map. The final averaged maps are reported in Figure 9, along with a set of representative snapshots for the metastable and other states. The DEK+1A FES shows a shallow minimum (z1) and two deep FE minima (z2 and z3). The z1 configuration is characterized by the two ions at the inner position, with their COMz at ∼8 Å and the lysine in proximity of the EEDD ring, nearly at the upper limit of the sampled space, rz(K-NH3+) ≈ 16 Å. The z2 and z3 states correspond to higher positions for the COMz of the ions and lower positions for K1422, i.e., COMz(Na1, Na2) ≈ 10 Å and rz(K-NH3+) ≈ 11 Å for z2 and COMz(Na1, Na2) ≈ 14 Å and rz(K-NH3+) ≈ 6 Å for z3. Hence, by considering the sequence of states in the order z3 to z1, we can describe the movement of two Na+ ions within the C/SF as follows: z3: the two ions are at their highest positions, while K1422 is at its lowest, and the sequence of positive charges from the extracellular side to the inner part of the filter is [Na+, Na+, K1422- N H3+]. z2: the ions move to a lower position (in particular the inner one), and K1422 switches to an intermediate one between them; the sequence of positive charges from the extracellular side to the inner region is now [Na+, K1422- N H3+, Na+]. z1: from the previous state, a less stable configuration is accessible, characterized by a slightly lower position of the ions and a considerably higher location of K1422; the sequence of positive charges is now [K1422- N H3+, Na+, Na+]. These multiple features are completely lost in the maps of the other three systems DEK0A, DEE–1A, and DEE0A. The DEK0A FES shows a single minimum, z4, defined by COMz(Na1, Na2) ≈ 8 Å and rz(K-NH2) ≈ 6 Å, i.e., with the ions and the uncharged lysine at lower positions. This configuration of the side chain is again similar to that observed in ref (70). When the DEKA filter is replaced by DEE–1A, the two ions are substantially trapped in the upper part of the SF (state z5, COMz(Na1, Na2) at ∼14 Å with one ion at the level of the EEDD motif and one ion in proximity of the DEEA ring), and their access to the inner region is precluded. In the case of the DEE0A system, the main metastable configuration (z6) corresponds to COMz(Na1, Na2) ≈ 9 Å and rz(E-COOH) ≈ 11 Å, with a less populated state with the uncharged glutamate at a lower position. This latter state is similar to z4, revealing again similarities among the two uncharged systems (DEK0A and DEE0A).

Voltage-gated sodium (Nav) channels play a pivotal role in the regulation of excitable cells, and mutations in genes encoding them have been linked to numerous disorders affecting nervous system function, heart rhythm, and muscle contraction. Therefore, understanding the structural and mechanistic features of the cationic fluxes through the Nav conductivity/selectivity filter (C/SF) is essential to accelerate the development of new efficient pharmacological treatments. In this context, a challenging but essential endeavor is to investigate the peculiar structural features of eukaryotic channels that are not obtainable from studying the high-resolution crystal structures of bacterial channels. In this study, we used a combination of molecular modeling, standard MD, and enhanced sampling to provide a description of ion permeation through the SF of human neuronal Nav channels. We combined results from MD simulations of homology-based models of Nav1.1, -1.2, and -1.6 and the recent cryo-EM Nav1.2 structure (PDB ID 6J8E) to describe the thermodynamics of Na+ residency within the C/SF. We ran simulations on multiple variants of the key SF DEKA ring by modifying the positively charged K (K1422 in Nav1.2), which is demonstrated to be essential for selectivity. Our simulations include a set of enhanced-sampling MD simulations employing the TAMD/OTFP method to reconstruct 2D FESs for single and double Na+ ions through the C/SF of the Nav1.2 structures. TAMD/OTFP has been previously assessed as an efficient tool to investigate systems with hidden FE barriers. FESs were obtained separately for two CV sets, the vertical position of one or two Na+ ions in the C/SF and that of the side chain of residue 1422. We also investigated the role of K1422 protonation and of charged and uncharged K1422E mutants. Moreover, the obtained FESs show that in the WT DEK+1A system there is a pathway through the C/SF with low FE values that the ions can follow to reach the inner portion of the filter, both as single- or double-ion configurations. We then described a putative mechanism of Na+ permeation in the case of the single-ion translocation using the string method. For all other systems (DEK0A, DEE0A, DEE–1A), the FESs display single isolated minima, mostly localized in the inner portion of the SF. Our simulations constitute a first attempt to investigate Na+ permeation in neuronal Nav channels employing mammalian cryo-EM structures. They clearly show the metastable sites for Na+ in the C/SF and the process by which they are accessed. Moreover, they reveal the role of the lysine residue in controlling Na+ translocation, in one- and two-ion occupancy states. This was investigated in previous work, albeit using a modified structure of a bacterial Nav channel where residues of the human Nav1.2 filter were grafted. The present approach, leveraging the TAMD/OTFP methodology, may be used in further works to extend the space of sampling of the simulations and, in addition, to explore the selectivity mechanisms of neuronal Nav channels. The use of cryo-EM structures in atomistic simulations requires caution, mostly because of the different resolution observed in the various protein domains. Even when good local resolution is available (as for the hNav1.4), the less resolved components of the structure can affect the fine balance between the residues lining the C/SF domain and the permeating ions. For this reason, during the simulations, it is important to monitor the structural stability and possible divergences from the initial configuration for the whole protein and for specific functionally relevant domains, as we reported here in detail. The definitive validation of cryo-EM structures and their characterization from a functional point of view remain long-term goals yet to be achieved. A significant step in this direction is represented by the recent boost of deep-learning-based approaches to predict protein structures, as it has been observed that predicted models of Nav channels are quite close to cryo-EM structures. However, limited differences have been pointed out in specific protein domains, and further experimental and computational studies are necessary to understand their effects. It is also important to comment on the use of homology/comparative modeling. In our case, the application of this strategy is supported by the high sequence identity between neuronal hNav channels and the hNav1.4 homologue (∼85% to ∼91% in the TM region, including the conserved EEDD and DEKA motifs in the C/SF), appropriate for obtaining high-quality homology models at atomic resolution. Such high sequence identity among the nine isoforms of hNavs has also functional consequences, and it is considered a cause for the limited selectivity of therapeutics that target these proteins, leading to various off-target side effects. As for the other protein domains, we employed a conservative approach and maintained the long extracellular loops that have been previously excluded in other studies. Our homology modeling relies on the use of state-of-the-art software for protein structure predictions, SWISS-MODEL, one of the most widely used and effective programs in the field. Moreover, the models were refined with FG-MD, a fragment-guided MD sampling developed by the I-TASSER community in order to improve the quality of multisubunit protein models. Finally, the recent cryo-EM Nav1.2 structure was also investigated, since it is the only experimental one available for the channel at the moment and is therefore of relevance to verify any conclusion drawn from homology models. Another key issue in atomistic simulations is the choice of the force field (FF). We employed here the most recent version of the CHARMM FF (CHARMM36m), which to the best of our knowledge is the most-used one by the ion channel community (including Navs) and membrane protein simulators. An interesting comparison between the effects of CHARMM and AMBER FFs was recently reported, showing marked differences in the relative probability of conductive/occluded configurations of the KcsA channel’s SF. These results leave open challenges in atomistic simulations, and the assessment of our results with other FFs will be the subject of future work.