Selective ion permeation involves complexation with carboxylates and lysine in a model human sodium channel

Bacterial and human voltage-gated sodium channels (Navs) exhibit similar cation selectivity, despite their distinct EEEE and DEKA selectivity filter signature sequences. Recent high-resolution structures for bacterial Navs have allowed us to learn about ion conduction mechanisms in these simpler homo-tetrameric channels, but our understanding of the function of their mammalian counterparts remains limited. To probe these conduction mechanisms, a model of the human Nav1.2 channel has been constructed by grafting residues of its selectivity filter and external vestibular region onto the bacterial NavRh channel with atomic-resolution structure. Multi-μs fully atomistic simulations capture long time-scale ion and protein movements associated with the permeation of Na+ and K+ ions, and their differences. We observe a Na+ ion knock-on conduction mechanism facilitated by low energy multi-carboxylate/multi-Na+ complexes, akin to the bacterial channels. These complexes involve both the DEKA and vestibular EEDD rings, acting to draw multiple Na+ into the selectivity filter and promote permeation. When the DEKA ring lysine is protonated, we observe that its ammonium group is actively participating in Na+ permeation, presuming the role of another ion. It participates in the formation of a stable complex involving carboxylates that collectively bind both Na+ and the Lys ammonium group in a high-field strength site, permitting pass-by translocation of Na+. In contrast, multiple K+ ion complexes with the DEKA and EEDD rings are disfavored by up to 8.3 kcal/mol, with the K+-lysine-carboxylate complex non-existent. As a result, lysine acts as an electrostatic plug that partially blocks the flow of K+ ions, which must instead wait for isomerization of lysine downward to clear the path for K+ passage. These distinct mechanisms give us insight into the nature of ion conduction and selectivity in human Nav channels, while uncovering high field strength carboxylate binding complexes that define the more general phenomenon of Na+-selective ion transport in nature.


Model Na v 1.2 Channel Creation
A model of the human Na v 1.2 channel was created by using the bacterial Na v Rh as a scaffold, for reasons outline in the main text Methods section. The structures of the two channels were aligned (see Fig.1) and the residues in and around the SF mutated into Na v Rh residues. The mutated region contained the key DEKA-ring as well as the vestibular EEDD-ring and surrounding residues believed to be important (see Methods). Models were created both with charged and neutral DEKA Lys, and proteins were then relaxed using Molecular Dynamics (MD) simulation in hydrated lipid bilayers, as detailed in the following section. Additional procedures not described in the main text Methods section are discussed below.
Rebuilding missing intracellular loops: Prior to sequence patching, modifications to the Na v Rh structure were required, owing to the absence of intracellular loops connecting S2 and S3 helices in the VSD. The Na v Rh structure was completed by rebuilding the missing loops (residue 70-75, Na v Rh numbering) in VSD II using Rosetta [1], chosen as it is the VSD with most residues resolved in the crystal structure of Na v Rh. This loop was then grafted on to the other VSDs, based on the assumption that these loops will be similar in all VSDs [2], minimizing the number of residues being rebuilt. The loop conformation chosen was representative of the best cluster in Rosetta and favors the positioning of aromatic residues F72 and F73 at the lipid interface, while charged residues D71 and K74 reside in the aqueous environment. Torsional angles of the connecting parts of the model were energy minimized without changing the side chain or backbone conformation of the rest of the protein. During loop optimization, the original X-ray structure rotamer orientations were used for residues not being rebuilt to keep the model as close as possible to the initial structure.
Maintaining key hydrogen bonds: There is a hydrogen bond present in both Na v Rh and Na v Ab structures between Thr178 and Trp182 that has been postulated to be important for keeping the shape of the SF in bacterial Na v [2,3]. These residues are conserved throughout the Na v family and all but one is present in human Na v 1.2 (T382-W386, C940-W943, T1420-W1424, T1712-W1716; Na v 1.2 numbering). To make sure that the mutations do not disturb the shape of the SF, this bond was therefore constrained throughout the mutation, minimization and equilibration of our model, and then released for unbiased production simulations.

Molecular Dynamics Simulation Details
The proteins were embedded in lipid bilayers of dipalmitoyl-phosphatidylcholine (DPPC), being the best characterized lipid for MD simulations [4], with explicit TIP3P water molecules [5] and 150 mM of NaCl or KCl solution. Systems were built and pre-equilibrated with CHARMM [6] and further equilibrated using NAMD [7] using the CHARMM36 [4] lipid and CHARMM22 protein and ion parameters [8] with CMAP corrections [9], to be consistent with previous simulations of the bacterial Na v Ab [10].
A number of attempts to adjust the extent of ion pairing between Na + and K + interacting with Glu/Asp carboxylate oxygen atoms have been made in the past by adjusting Leonard-Jones (LJ) parameters with a non-bonded fix (NBFIX) according to hydration free energy and osmotic pressure data [11,12,13]. These investigations had an overall similar conclusions on ion selectivity when using modified and unmodified parameters for LeuT [11], a Na + /Ca 2+ exchanger [12], and an acid sensing ion channel [14]. Furthermore, previous simulations of ion pairing between Na + and acetate have shown reasonable agreement between standard LJ parameters and quantum mechanical ab initio MD Umbrella Sampling calculations, and suggest that addition of NBFIX adjustment makes the ion-carboxylate pairing far too repulsive, despite improving the description of osmotic pressure [14]. The binding affinities and dissociation constants (K D ) using either the unmodified Lorentz-Berthelot mixing rule or NBFIX modifications showed general agreement with experiment [14]. We have previously shown, as expected, that changing to a more repulsive NBFIX interaction between ions and carboxylates leads to a reduction in ion-carboxylate pairing, but that even then Na + maintains increased multi-ion complexes compared to K + [14]. Therefore, we have used standard parameters for the ion-carboxylate interaction in this study, and refer to Ref. [14] for more information on the effects of NBFIX on ioncarboxylate complex formation.
Adjustments to the LJ NBFIX for ions interacting with backbone carbonyl oxygen atoms ( ! =−0.102 and !"# =3.64 Å for K + -O and ! =−0.0750 and !"# =3.30 Å for Na + -O) have previously been made to yield free energies of ion solvation (of -80 kcal/mol and -101 kcal/mol for K + and Na + ) in the protein backbone mimetic, Nmethyl-acetamide that match free energies of ion hydration [15]. These parameters have been previously used for narrow K + channels [16] and also adopted for these simulations to better describe Na + and K + interacting with backbone carbonyl oxygen atoms (during Anton production runs, described below, these were only applied to SF residues 178-181 due to a limitation in number of NBFIX modifications permitted).
During equilibration simulations using NAMD 2.9 [17], the NPT ensemble was imposed using a Langevin piston to maintain constant pressure of 1atm [18,19], and a Nosé-Hoover thermostat to maintain constant temperature of 323K [20,21]. The RATTLE algorithm [22] was used to maintain the bonds to hydrogen atoms. Electrostatic interactions were calculated for the tetragonal periodic box using Particle Mesh Ewald [23] with grid spacing of 1.5 Å and 6 th order B-spline for mesh interpolation. The non-bonded pair list was built to 16 Å, with a real space cut off at 12 Å and force switching from 10 Å. The protein was constrained while the membrane was energy minimized for 1000 steps and then equilibrated by initially constraining all heavy atoms with a 10 kcal/mol/Å 2 harmonic restraint, which was then slowly released over 2.5 ns. Residue 182 in subunit II was removed and residue 181 and 183 pulled together in steps of 0.1 Å with a force constant of 100 kcal/mol/Å 2 , decreased over 20 ns, while the other subunits and residues outside of the SF were constrained during this adjustment, with protein heavy atoms constrained using 10 kcal/mol/Å 2 harmonic constraints. Restraints with force constant 10 kcal/mol/Å 2 on the hydrogen bonds between T382-W386, C940-W943, T1420-W1424 and T1712-W1716 were added, and the protein was equilibrated for an additional 100 ns. Unbiased production simulations were then carried out on the purpose-built supercomputer Anton [24,25] for 4 for Na v Rh/Na v 1.2 with charged Lys both for NaCl and KCl, and 2 for neutral Lys (less time required due to better sampling of ion movements), totaling 12 . The NPT ensemble was maintained on Anton with a Nose-Hoover thermostat [20,21] and a Martyna-Tobias-Klein barostat [26] at 323 K with a 2 fs time step. A cut-off of 12.72 Å was used and long-range electrostatic forces evaluated with Gaussian split Ewald method every 6 fs [27,28,29].

Simulation Analysis
To determine the number of ions involved in ion conduction, the occupancy in the pore was calculated based on the mean and standard deviation of the number of ions within r < 10 Å and −15 < z − z !"# < 15Å, where z ref is the z position of the center of mass (COM) of the backbone atoms of residues 178-182 (incorporating the DEKA as well as two additional residues above and two residues below); chosen to encompass entire Na v Rh/Na v 1.2 pore, including accessible region of the intracellular gate, central cavity, SF and the outer vestibule. This outer cut off of z=+15 Å is consistent with our analysis of the vestibular ring structure, such that this range always includes the EEDD side chains. Compared to previous studies with Na v Ab [10], the outer cut off was extended by 2.5 Å to allow for this broader binding region in the outer vestibule. Permeation events were counted every time an ion crossed the SF and passed the Lys (inward moving ion starting at > 5 Å and reaching < !"# and < −2 Å; reversed for an outward permeating ion).
The free energy map, or potential of mean force (PMF), for ion movement through the filter was calculated for each system using trajectory data with specific ion occupancies, to ensure well-defined free energy maps. The relevant occupancies for permeation were found to be either 2 or 3 ions for the channel with neutral Lys, and 1 or 2 ions for the channel with charged Lys. The PMFs were calculated from unbiased simulations as where ρ is the unbiased probability distribution as a function of reaction coordinate(s) ! , being the vertical positions of one or more ions (z 1 , z 2 , or z 3 ; defined with z 1 being the lower-most ion), or their centroids (e.g. !" = ( ! + ! )/2 ), each relative to !"# , and where C is a constant. Errors for the maps were calculated by dividing the well-sampled part of the trajectory (the last 2 s of the trajectory for NaCl and the last 2.2 s trajectory for KCl) into two parts to determine the deviation between the maps. The total error for the maps was calculated as the root mean square, including all points with energy < 3 kcal/mol (to avoid poorly sampled regions at the periphery).
Structural changes, reported in S4 Fig, were investigated by root mean square deviation (RMSD) of the entire mutated region (residues 174-184). The RMSD of the mutated region was calculated to determine the stability of the model. S4 Fig shows  how, although the model goes through initial structural changes, the RMSD becomes stable and changes subside. Experiments have shown that the SF of the mammalian Na v is highly flexible [30], placing some limitation on direct structural comparisons between model and published eukaryotic structures. Therefore, the Root Mean Square Fluctuation (RMSF) was calculated for the backbone of each residue in the SF and vestibule (spanning only the SF residues 178-184) to evaluate the flexibility of the SF region only, reported in S6 Fig. The RMSF was measured after orientation with respect to the backbone of residues 178-184 for either all subunits together, or each individual subunit, to investigate both relative and internal fluctuation, respectively. For comparative purposes, the same calculation was done for 2.5 s long simulations with the bacterial channel Na v Ab. The Na v Rh/Na v 1.2 model was also compared structurally to the eukaryotic Na v PaS and EeNa v 1.4 by performing RMSD calculations of the equilibrated Na v Rh/Na v 1.2 model to the cryo-EM structure of Na v PaS and EeNa v 1.4, using the backbone of each residue in the SF and vestibule (residue 178-184). The stability of the H-bonds between T382-W386, C940-W943, T1420-W1424 and T1712-W1716 (Na v 1.2 numbering, corresponding to Thr178-Trp182 in Na v Rh) were evaluated by calculating the donor-acceptor distance over the course of the simulation.
To determine the important states involved in ion permeation, clustering analysis was performed. The frames were broken down according to ion occupancies and the representative occupations further analyzed. The radial distribution function for ionion and ion-carboxylate pair interactions was evaluated and a cutoff based on the first minimum was then determined (Na + -Na + :r min =4.7 Å, K + -K + : r min =5.5 Å, Na +carboxylate:r min =3.8 Å, K + -carboxylate:r min =4.2 Å) (S3 Fig). This cutoff was used to assign configurations to clusters based on the number of DEKA and EEDD carboxylates within the ion-carboxylate cutoff and, in the case of multi-ion occupancies, whether or not the ions were within the ion-ion cutoff distance. The clusters were further broken down according to which residues the ions were bound to. The most stable frame of the most representative clusters was identified by choosing the configuration with highest probability of occurrence within each cluster.
The most important states involved in the permeation of Na + and K + ions across the SF were identified from the PMFs (Fig.3&4) and cluster analysis (S1 Table). These states include complexes involving one or more carboxylate groups, Na + or K + ions and/or the signature Lys ammonium group. To evaluate the relative binding affinities of these states, Free Energy Perturbation (FEP) [31] calculations were performed. The starting frame was chosen as the most stable frame of the cluster, as described above. The relative binding affinity, !! ! →! ! , was calculated from The NAMD FEP module, using dual topology, was used to determine !! ! →! ! :!"#$ and !! ! →! ! :!"#$ by alchemically transforming Na + to K + and K + to Na + using a coupling parameter, . The coupling parameter was divided into 20 windows that were initiated by consecutively increasing in steps of 0.05 for 0.2 ns each. Each window was then equilibrated and simulated concurrently for 2.4 ns giving a total of 48 ns simulation per calculation. The windows were combined and the free energy difference calculated with WHAM [32]. Both forward and backward transformations were performed (Na + →K + and K + →Na + ), and for the case of a site bound by 2 ions, each individual ion was transformed separately and the final transformation calculated by the mean of the two transformations (1 st ion: mean of Na + Na + →K + Na + and Na + Na + → Na + K + , or mean of K + K + →K + Na + and K + K + → Na + K + ; 2 nd ion: mean of K + Na + →K + K + and Na + K + → K + K + ; or mean of K + Na + →Na + Na + and Na + K + → Na + Na + ). The endpoint histograms for the adjacent for the forward and backward transformations were combined with the VMD plugin ParseFEP [33] using Bennett Acceptance Ratio (BAR) [34].
During these FEP calculations, ions were constrained with a 20 kcal/mol/Å 2 halfharmonic flat-bottomed constraint with outer radius of 2.75 Å from the COM of the multiple carboxylate carbon atoms for the ion/ion+Lys-carboxylate clusters. The results were adjusted for the constraint force by subtracting the Boltzmann average of the difference in constraint energy for = 0 and = 1. Most of the time the constraints were not interfering and thus only a small correction (no greater than 0.3 kcal/mol; computed as the difference in − ! exp (−Δ / ! ) (! ! (!)) between = 0 and 1, where U 0 ( ) is the unperturbed system potential and Δ is the constraint energy) was subtracted from the results. For the case of ion+Lys in the lower SF, a constraint was applied to the ion and the ammonium group of the Lys, keeping them at that location (-4.0<z<-2.5 Å), using a flat-bottomed harmonic constraint with force constant 20 kcal/mol/Å 2 . The correction for the constraint was found to be 0.34±0.12 kcal/mol. During the Na + → K + and K + → Na + transformations, a harmonic constraint of 100 kcal/mol/Å 2 was applied between the vanishing and appearing ions to maintain overlap. A convergence criterion of within 1 kcal/mol was used. The first 1.2 ns of each window was discarded for the single-ion/multi-carboxylate complex. For the 2-ion complex, the first 0.6 ns for the 1 st and 1.4 ns for the 2 nd ion transformation were removed from each window. For the ion+Lys complex, the first 1.0 ns was discarded, while for the ion adjacent to Lys in the lower SF the first 1.2 ns was discarded for each window.
To examine the electrostatic influence of SF signature sequences and vestibular charge rings on permeation, Poisson's equation was solved for long ranged electrostatic interactions using a smooth particle mesh Ewald method, using VMD plugin PMEPot [35], yielding electrostatic potential maps for the Na v Rh/Na v 1.2 and Na v Ab channels. A grid size of 1 Å was used and the final map in the SF was averaged over every 100th frame. The ion density was evaluated with the VMD VolMap plugin [36].