Determination of the μ-Conotoxin PIIIA Specificity Against Voltage-Gated Sodium Channels from Binding Energy Calculations

Voltage-gated sodium (NaV) channels generate and propagate action potentials in excitable cells, and several NaV subtypes have become important targets for pain management. The μ-conotoxins inhibit subtypes of the NaV with varied specificity but often lack of specificity to interested subtypes. Engineering the selectivity of the μ-conotoxins presents considerable complexity and challenge, as it involves the optimization of their binding affinities to multiple highly conserved NaV subtypes. In this study, a model of NaV1.4 bound with μ-conotoxin PIIIA complex was constructed using homology modeling, docking, molecular dynamic simulations and binding energy calculations. The accuracy of this model was confirmed based on the experimental mutagenesis data. The complex models of PIIIA bound with varied subtypes of NaV1.x (x = 1, 2, 3, 5, 6, 7, 8, or 9) were built using NaV1.4/PIIIA complex as a template, and refined using molecular dynamic simulations. The binding affinities of PIIIA to varied subtypes of NaV1.x (x = 1 to 9) were calculated using the Molecular Mechanics Generalized Born/Surface Area (MMGB/SA) and umbrella sampling, and were compared with the experimental values. The binding affinities calculated using MMGB/SA and umbrella sampling are correlated with the experimental values, with the former and the latter giving correlation coefficient of 0.41 (R2) and 0.68 (R2), respectively. Binding energy decomposition suggests that conserved and nonconserved residues among varied NaV subtypes have a synergistic effect on the selectivity of PIIIA.


Introduction
The µ-conotoxins are short snail peptide toxins specifically targeting subtypes of Na V and are considered as important drug leads. The first step for engineering the specificity of µ-conotoxins is understanding of the molecular interaction mechanism between the µ-conotoxins and varied subtypes of the Na V . Accurate determination of the molecular determinants that confer the specificity of the µ-conotoxin PIIIA to different Na V subtypes is essential for further engineering the specificity of the µ-conotoxins for therapeutic purposes.  (1)(2)(3)(4) were colored in green. The segments 5 and 6, as well as the sequence between TM5 and TM6 (P1 loop) (Display in the B), that enclose to form the pore domain of the NaV were colored in cyan, orange, and magenta, respectively; (B) The cryo-EM structure of a eukaryotic NaV channel (PDB code: 5X0M). The structure contains four voltage-sensing domains and one pore domain. TMD represents transmembrane domain, while TM5 and TM6 represent the segments 5 and 6 of transmembrane domain; (C) Crystal structure of the prokaryotic sodium channel from Magnetococcus marinus (NaVMs) (PDB code: 4OXS). The structure of the NaVMs was colored according to the topology of the eukaryotic NaV channel.  (1)(2)(3)(4) were colored in green. The segments 5 and 6, as well as the sequence between TM5 and TM6 (P1 loop) (Display in the B), that enclose to form the pore domain of the NaV were colored in cyan, orange, and magenta, respectively; (B) The cryo-EM structure of a eukaryotic Na V channel (PDB code: 5X0M). The structure contains four voltage-sensing domains and one pore domain. TMD represents transmembrane domain, while TM5 and TM6 represent the segments 5 and 6 of transmembrane domain; (C) Crystal structure of the prokaryotic sodium channel from Magnetococcus marinus (NaVMs) (PDB code: 4OXS). The structure of the Na V Ms was colored according to the topology of the eukaryotic Na V channel.

Conotoxins
Conotoxins are short disulfide-rich peptide toxins comprising 12-30 amino acid residues and possess varied selectivity on ion channels [9]. µ-Conotoxins are selective and active blockers of the voltage-gated sodium channels. The µ-conotoxin PIIIA was isolated from Conus purpurascens, an Eastern Pacific fish hunting species of cone snails [10]. The PIIIA possesses varied potency towards TTX sensitive Na V channels [11,12], and therefore the PIIIA serves as a valuable probe to distinguish Na V channel subtypes when used in conjunction with other selective toxins [11,13] (Figure 2).

Conotoxins
Conotoxins are short disulfide-rich peptide toxins comprising 12-30 amino acid residues and possess varied selectivity on ion channels [9]. μ-Conotoxins are selective and active blockers of the voltage-gated sodium channels. The μ-conotoxin PIIIA was isolated from Conus purpurascens, an Eastern Pacific fish hunting species of cone snails [10]. The PIIIA possesses varied potency towards TTX sensitive NaV channels [11,12], and therefore the PIIIA serves as a valuable probe to distinguish NaV channel subtypes when used in conjunction with other selective toxins [11,13] (Figure 2). Without the high-resolution crystal structure of the eukaryotic voltage-gated sodium channel, crystal structures of the prokaryotic sodium channel NaVAb have been widely used as templates to model the interactions between μ-conotoxins and NaV1.4 [14,15]. These models were used to predict the binding affinities of μ-conotoxins against NaV1.4, and explained their interaction mechanism [14,15]. However, the molecular determinants that confer the specificity of the μ-conotoxins to varied Without the high-resolution crystal structure of the eukaryotic voltage-gated sodium channel, crystal structures of the prokaryotic sodium channel Na V Ab have been widely used as templates to model the interactions between µ-conotoxins and Na V 1.4 [14,15]. These models were used to predict the binding affinities of µ-conotoxins against Na V 1.4, and explained their interaction mechanism [14,15]. However, the molecular determinants that confer the specificity of the µ-conotoxins to varied Na V subtypes are still elusive, and no effective computational methodologies are available to predict the specificity of the µ-conotoxins.
In this study, the molecular determinants that confer the specificity of PIIIA to varied subtypes (hNa V 1.1-1.9) of the Na V were determined through extensive molecular dynamics simulations and binding energy calculations. And a computational method that could be used to predict the specificity of the µ-conotoxin PIIIA against different Na V subtypes was proposed.

Flowchart for Specificity Prediction
Despite millions of years of divergent evolution between the Na V Ms and eukaryotic sodium channel, they both still possess similar pharmacological profiles to the channel blockers as shown by the functional electrophysiological studies [16]. Indeed, previous computational modeling studies suggested that crystal structure of NaChBac (Na(+)-selective channel of bacteria) could be used as template to model the pore domain of the human or rat Na V 1.4 due to the relatively high sequence identity and structural similarity at the SF domain of the channel. In this study, the crystal structure of the Na V Ms in open (PDB code 4CBC) was used as template to model the pore domain of the human Na V 1.4. After building the homology model of Na V 1.4, the NMR structure of the PIIIA (PDB code 1R9I) was docked into the extracellular pore domain of the Na V 1.4 model [17]. Ten conformations of the PIIIA with minimum docking score were maintained and subjected to further molecular dynamics (MD) refinement and binding energy calculations using MMGB/SA method. The Molecular Mechanics/Poisson Boltzmann Surface Area (MM/PBSA) and Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) approaches are more computationally efficient which estimate protein-protein binding free energies. Conformations of the PIIIA at Na V 1.4 were ranked based on the binding energies calculated using MMGB/SA. And relatively higher energy conformations of the PIIIA at Na V 1.4 were discarded, whereas the low energy conformations were maintained and further validated based on the known mutagenesis data [10,11]. Once the determined binding mode of PIIIA prokaryotic Na V 1.4 was validated by the experimental data, it could be used as template to model interactions of PIIIA with other Na V subtypes.
The binding modes of PIIIA against Na V 1.
x (x = 1, 2, 3, 5, 6,7,8,9) were determined based on homology modeling using PIIIA/Na V 1.4 as the template. The produced models were further refined using MD simulations. The binding energies of PIIIA against Na V 1.x (x = 1, 2, 3,5,6,7,8,9) were calculated using MMGB/SA method and umbrella sampling. The final step is the analysis of the correlation between the calculated and experimental derived binding affinities. The flowchart for PIIIA specificity prediction was summarized in Figure 3.

Homology Modeling of the Na V 1.4
Models of the Na V 1.4 were built using the crystal structure of prokaryotic Na V Ms channel (PDB code: 4P9O) as template in MODELLER (version 9v12) [18,19], as described previously [20]. The sequence of Na V 1.4 was retrieved from the UniProt database [21], and was divided into four homologue sequence fragments based on the sequence of prokaryotic Na V Ms. The sequence alignment between the four homologue sequence fragments of Na V 1.4 and sequence of the prokaryotic Na V Ms ( Figure S1) was used to build 100 homology models of the Na V 1.4, and the model with the lowest dope score [22] was selected for PIIIA docking.

PIIIA Docking
Docking the NMR structure of PIIIA (PDB Code: 1R9I) to NaV1.4 model was performed using AutoDock 4.2 [23]. Gasteiger charges were used, and nonpolar hydrogens of the macromolecule and ligand were merged. A grid box with dimensions of 60 Å × 60 Å × 60 Å and a grid spacing of 0.375 Å was set up and centered on the geometry center of the "P1 loop" of the four domains. Docking was performed using a Lamarckian genetic algorithm (LGA), with the receptor treated as rigid, the docking results were analyzed by AutoDock Tools. The produced conformations with top 10 docking score were selected for the MD refinement.

Restraint Molecular Dynamics Simulations
The selected 10 NaV1.4/PIIIA complexes were subjected to MD simulations for refinement of the PIIIA binding mode in AMBER 16 [24]. The MD simulations are performed using the ff14SB force field for the proteins and the peptide [25]. Parameters for the pyroglutamic acid (PCA) residue were prepared using the antechamber module in AMBER 16. Atom partial charges for PCA (pyroglutamic acid) were produced using the R.E.D Tools [26]. The PCA was incorporated into the PIIIA peptide as an unnatural amino acid. After building the whole system in LEaP that is a module of AMBER16, Figure 3. Flowchart for prediction of the specificity of the µ-conotoxin PIIIA to varied Na V subtypes.

PIIIA Docking
Docking the NMR structure of PIIIA (PDB Code: 1R9I) to Na V 1.4 model was performed using AutoDock 4.2 [23]. Gasteiger charges were used, and nonpolar hydrogens of the macromolecule and ligand were merged. A grid box with dimensions of 60 Å × 60 Å × 60 Å and a grid spacing of 0.375 Å was set up and centered on the geometry center of the "P1 loop" of the four domains. Docking was performed using a Lamarckian genetic algorithm (LGA), with the receptor treated as rigid, the docking results were analyzed by AutoDock Tools. The produced conformations with top 10 docking score were selected for the MD refinement.

Restraint Molecular Dynamics Simulations
The selected 10 Na V 1.4/PIIIA complexes were subjected to MD simulations for refinement of the PIIIA binding mode in AMBER 16 [24]. The MD simulations are performed using the ff14SB force field for the proteins and the peptide [25]. Parameters for the pyroglutamic acid (PCA) residue were prepared using the antechamber module in AMBER 16. Atom partial charges for PCA (pyroglutamic acid) were produced using the R.E.D Tools [26]. The PCA was incorporated into the PIIIA peptide as an unnatural amino acid. After building the whole system in LEaP that is a module of AMBER16, minimization was performed to remove the van der Waals clashes between PIIIA and Na V 1.4. The 2000 steps of steepest descent minimization and 3000 conjugate gradient minimum were performed with the solute restrained using a harmonic force potential with a spring constant of 100 kcal/mol −1 /Å −2 . After the first round of minimization, the restraints were withdrawn and the whole system was minimized using the same parameters as above. MD simulations were carried out after minimization. First, the whole system was gradually heated from 50 K to 300 K at constant volume and temperature ensemble for 100 ps with the solute restrained with a harmonic force of 5 kcal/mol −1 /Å −2 . The simulations were then switched to constant pressure and temperature ensemble by maintaining the harmonic restraints in 100 ps. For the production run, the restraints on PIIIA as well as the P1 loop were withdrawn and 20 ns MD simulations were performed with the temperature and pressure maintained at 300 K and 1 bar, respectively. All bonds involving hydrogen atoms were constrained with the SHAKE algorithm with a time step of 2 fs [27]. The particle-mesh Ewald (PME) method was used to model long-range electrostatic interactions [28]. The following MD simulations on PIIIA/ Na V 1.x (x = 1, 2, 3, 5, 6,7,8,9) used the same parameters setting up as PIIIA/Na V 1.4 system.

Molecular Dynamics Simulation of the Na V 1.4 in Membrane
For evaluation of the membrane influences on the conformation of the Na V 1.4 and the bound PIIIA, MD simulations were performed on Na V 1.4 in apo and in complex with PIIIA. The Na V 1.4 in apo or bound with PIIIA were inserted a bilayer containing a 2:2:1 mixture of POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine):POPE (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine):Cholesterol as suggested in previous study [29] with dimension of 80 × 80 × 80 Å 3 and the system solvated with TIP3P [30] water molecules and Na + and Cl − ions such that the system was neutral with an overall concentration of 0.15 M in CHARMM-GUI that is a webserver for preparing MD simulation system (http://www.charmm-gui.org). The temperature of the system was gradually increased to 310 K and was equilibrated for 500 ps in NVT (constant-volume, constant-temperature) and NPT (constant-pressure, constant-temperature) ensembles, respectively, with the protein and lipids restrained with 10 kcal/mol/Å 2 force. Langevin thermostat is used for the initial heating. For the second phase heating, an anisotropic Berendsen weak-coupling barostat is used to equilibrate the pressure in addition to the use of the Langevin thermostat to equilibrate the temperature. Then, the restraint of the membrane was withdrawn and the whole system was simulated for 20 ns in NPT for properly equilibrating the membrane system. And the restraint on the protein was gradually withdrawn in 10 steps in 5 ns MD simulations. Afterwards, 150 ns production run was performed on the Na V 1.4 in apo and in complex with PIIIA, respectively. In production run, the temperature was controlled using the Langevin thermostat while pressure was controlled using the anisotropic Berendsen barostat. All simulations were performed using the Lipid14 force field [31] for the lipid, AMBER14SB force field for the protein [31]. The MD simulations used a time step of 2 fs and, all bonds involving hydrogen atoms were maintained to their standard length using the SHAKE algorithm [32]. Particle Mesh Ewald (PME) [33] was used with a cutoff of 10 Å for non-bonded atoms interactions and neighbor lists were updated every 10 steps.

Binding Energy Calculations Using MMGB/SA
Molecular mechanics generalized Born surface area (MMGB/SA) [34] were applied to calculate the binding affinities of PIIIA against the Na V channels. The energies were averaged on the 50 frames extracted from the last 10 ns MD simulations. The parameters are as previously reported [20]. Briefly, the internal dielectric and external dielectric constants were set to 2.0 and 80.0, respectively. A probe radius of 1.4 Å, a grid spacing of 0.5 Å and ionic strength 0.1 mol/L were set up for the calculations.

Decomposition of the Binding Energies
The binding energies calculated using MMGB/SA were decomposed into energy contributions from each residue of the binding site using the MMPB/SA.py script in AMBER16. Parameters setting up was the same as above.

Binding Energies Calculations Using Umbrella Sampling
Umbrella sampling is one biased MD method that provides free energy along a reaction coordinate [35]. In umbrella sampling, a bias potential is applied to drive a system from one thermodynamic state to another along a reaction coordinate. Here, the distance between the center of mass of the PIIIA and the center of mass of the five symmetric CA atoms in the transmembrane domain were selected as the reaction coordinate, and along which a harmonic potential was added (Equation (1)).
where the k is 100 kcal/mol/Å 2 and x 0 is the initial distance. When the reaction coordinate (x − x 0 ) reached 15 Å, the PIIIA was considered completely unbound from the channel. The reaction coordinate was divided into 30 windows, and in each window, 2 ns MD simulations were performed. The windows were then combined by the weighted histogram analysis method (WHAM) [36]. Conformations in the last 1.5 ns were used to derive the potential of mean force along the coordinate using WHAM. The binding/unbinding free energy of PIIIA was calculated using the PMF values at the end of unbinding.  Figure S2A,B). The binding energies predicted using Autodock are considerably higher than the experimental values. The unfavorable binding energies given by the docking technique are originated from the poor accuracy when dealing with the flexibility of the protein and the peptide by AutoDock. During docking, only side chains of the PIIIA were considered as flexible, whereas side chains of the Na V 1.4 were considered as rigid. As a consequence, the predicted binding affinities were underestimated using docking. Thus, it is prerequisite to refine the binding mode of PIIIA at Na V 1.4 using MD, and re-rank the binding affinities of PIIIA at Na V 1.4 using more rigorous energy estimation method, like MMGB/SA or MMPB/SA.

Results and Discussion
Here, the MMGB/SA method was used to calculate the binding affinities of the 10 conformations of PIIIA at Na V 1.4 due to its slightly better performance than MMPB/SA method to rank the binding affinities of the α-conotoxin analogues in a previous study [20]. The ranking order for the 10 conformations significantly changed after re-ranking based on the MMGB/SA method (Table S1), and their backbone orientation at the binding site of Na V 1.4 in some cases significantly varied ( Figure S3). Conformation for the side chains of PIIIA varied significantly after MD simulations in most cases (not shown).
Conformations of PIIIA with optimum binding energy were selected for accuracy validation by comparison with the mutagenesis data [37]. We found that the minimum binding energy conformation of PIIIA was the most consistent with the mutagenesis data, and its binding mode to Na V 1.4 was shown in Figure 3. The PIIIA was positioned at the pore surface of the ECD with the axis of the helix between residue 6 and 15 tilted to the pore, resulting in the N-termini outward tilted from the binding site forming no direct interactions with residues of the binding site. The predicted binding mode of PIIIA is consistent with the previously determined binding mode of PIIIA at the Na V 1.4 using restraint docking method [14]. The PIIIA is rich with positively charged residues, whereas the binding site of Na V 1.4 possesses an abundance of negatively charged residues ( Figure S4). Thus, the net charge of the PIIIA is overall complementary to the Na V 1.4. Indeed, there are numerous salt-bridges formed at the interface between PIIIA and Na V 1.4 ( Figure 4). The R14 and K17 of PIIIA form salt-bridges with E755 and D1241, two of the filter residues at the pore domain, respectively. The R12 of PIIIA forms salt-bridges with E758 that is near the pore filter, while the R20 forms salt-bridges with D1541 at the periphery of the binding site. Besides, the hydroxyl group of O18 forms a pair of hydrogen bonds with N1536. Thus, our predicted binding mode of PIIIA against the Na V 1.4 suggests that the salt-bridges as well as the hydrogen bonds are the dominant binding determinants for PIIIA binding to Na V 1.4. Indeed, the previous mutagenesis data also supports that R14, R12, K17, and R20 contribute most to the binding affinity of PIIIA to Na V 1.4 [37]. More detailed pairwise interactions between PIIIA and Na V 1.4 were listed in Table 1. O18 forms a pair of hydrogen bonds with N1536. Thus, our predicted binding mode of PIIIA against the NaV1.4 suggests that the salt-bridges as well as the hydrogen bonds are the dominant binding determinants for PIIIA binding to NaV1.4. Indeed, the previous mutagenesis data also supports that R14, R12, K17, and R20 contribute most to the binding affinity of PIIIA to NaV1.4 [37]. More detailed pairwise interactions between PIIIA and NaV1.4 were listed in Table 1.

Influences of the Membrane to the Conformation of the NaV1.4 in apo and NaV1.4 in Complex with PIIIA
For consideration of the influences of the membrane to the conformation of NaV1.4 and PIIIA/NaV1.4 complex model, nonrestraint molecular dynamic simulations were performed on the NaV1.4 model in apo or in complex with PIIIA embedded in the membrane. MD simulations results suggest that conformation of the NaV1.4 model in apo and in complex with PIIIA is stable in 100 ns MD. The RMSD (root-mean-square deviation) value of the NaV1.4 model in complex with PIIIA is lower than that in apo ( Figure S5A). Obviously, binding PIIIA to NaV1.4 has a stabilization effect on this model. After 150 ns of non-restraint MD the binding mode of PIIIA at NaV1.4 is essentially the same to that determined using restraint MD without adding the membrane ( Figure S5B). Thus, it is applicable to use more efficient restraint MD for the refinement of the binding mode of PIIIA at NaV1.4 without incorporation of the membrane.

Influences of the Membrane to the Conformation of the Na V 1.4 in apo and Na V 1.4 in Complex with PIIIA
For consideration of the influences of the membrane to the conformation of Na V 1.4 and PIIIA/Na V 1.4 complex model, nonrestraint molecular dynamic simulations were performed on the Na V 1.4 model in apo or in complex with PIIIA embedded in the membrane. MD simulations results suggest that conformation of the Na V 1.4 model in apo and in complex with PIIIA is stable in 100 ns MD. The RMSD (root-mean-square deviation) value of the Na V 1.4 model in complex with PIIIA is lower than that in apo ( Figure S5A). Obviously, binding PIIIA to Na V 1.4 has a stabilization effect on this model. After 150 ns of non-restraint MD the binding mode of PIIIA at Na V 1.4 is essentially the same to that determined using restraint MD without adding the membrane ( Figure S5B). Thus, it is applicable to use more efficient restraint MD for the refinement of the binding mode of PIIIA at Na V 1.4 without incorporation of the membrane.

Binding Modes of PIIIA at Other Subtypes of Na V
PIIIA possesses broad inhibitory spectra on Na V , from Na V 1.1 to Na V 1.8 [10][11][12][38][39][40]. Residues at the binding site of the α subunit of Na V is highly conserved, and it might be applicable to homology model the binding modes of PIIIA to Na V 1.x (x = 1, 2, 3, 5, 6, 7, 8, 9) by using PIIIA/Na V 1.4 as the template. In this study, the PIIIA was directly modelled into the binding site of Na V 1.x (x = 1, 2, 3, 5, 6, 7, 8, 9) using PIIIA/Na V 1.4 as a template instead of using the docking method in considering of its incapacity for dealing with the flexibility of the protein and the peptide and its inaccuracy for prediction the binding affinity of the peptide ligands. Then the binding modes of PIIIA to Na V 1.x (x = 1, 2, 3, 5, 6,7,8,9) were refined using MD simulations.
Binding modes of PIIIA to Na V 1.

Prediction of the Specificity of PIIIA to Varied NaV Subtypes
Prediction of the specificity of PIIIA to different Na V subtypes presents a challenge, as it involves in differentiation of the binding affinity of PIIIA to several highly conserved Na V subtypes. In this study, we used two computational methods, MMGB/SA and umbrella sampling to predict the binding affinities of PIIIA to Na V subtypes from Na V 1.1 to Na V 1.9.
The MMGB/SA was used to predict the binding affinities of PIIIA to varied Na V subtypes due to its better performance than MMPB/SA at ranking the binding affinities of α-conotoxin analogues in our previous study [21]. For sufficiently sampling the conformation of PIIIA at varied Na V subtypes, the top three DOPE score PIIIA/ Na V 1.x (x = 1, 2, 3, 5, 6, 7, 8, 9) models were selected for MD simulations, and 20 ns MD simulations were performed on each of the PIIIA/Na V 1.4 systems. Frames from the last 10 ns MD were used to calculate the binding energy using MMGB/SA method. The calculated binding affinities were averaged on the three MD trajectories for each of PIIIA/Na V 1.x (x = 1, 2, 3,4,5,6,7,8,9) systems. In a previous study, we found that incorporation of the entropy component to the G-free energy gave slightly worse ranking of the binding affinity of the α-conotoxin analogues despite of the significantly increased consumption of the computational resources [21]. Thus, only the enthalpy rather than the entropy was calculated and used for ranking the binding affinity of PIIIA to varied Na V subtypes. As shown by Table 2, the predicted binding energies without entropy using MMGB/SA are much lower than the experimental values, whereas they are modestly correlated with the experimental values giving a correlation coefficient of 0.41 (R 2 ) ( Figure 6). Thus, MD simulations coupling of MMGB/SA energy calculation could give modestly accurate predictions of the specificity of PIIIA to varied Na V subtypes.
The umbrella sampling is more efficient than MD simulations for conformation sampling, and it has been frequently applied to predict the binding free energies of the conotoxins to their target proteins [41][42][43][44][45][46]. In this study, umbrella sampling was used to predict the binding affinities of PIIIA against varied Na V subtypes. The above MD refined the top three models of the PIIIA/Na V systems, and were used as the initial for umbrella sampling simulations. The predicted binding free energies of PIIIA against varied Na V subtypes were derived from the calculation of potential of mean force (PMF) and were shown in Table 2. The absolute values of the predicted binding free energies using umbrella sampling are quite similar to the experimental values, and are high correlated to the experimental values with a correlation efficient of 0.68 (R 2 ). Thus, umbrella sampling could give a highly accurate ranking of the inhibition activity of PIIIA to different Na V subtypes. Umbrella sampling gave a better prediction of the specificity of PIIIA to varied Na V subtypes due to its more efficient conformation sampling. The binding affinities of PIIIA to varied Na V subtypes were calculated using MMGB/SA method. Twenty MD simulations were performed on the top 3 PIIIA/Na V models from homology modeling. The MMGB/SA calculation was performed on the last 10 ns of the five MD trajectories; b The binding affinities of PIIIA were calculated using umbrella sampling (US) on the top 3 PIIIA/Na V models; c Binding affinities of PIIIA were derived from the IC 50 of PIIIA to varied Na V subtypes from Na V 1.1 to Na V 1.9 using ∆G = −RTln(IC 50 ). Binding affinities of PIIIA were derived from the IC50 of PIIIA to varied NaV subtypes from NaV1.1 to NaV1.9 using ∆G = −RTln(IC50).

Binding Energy Decomposition Using MMGB/SA
The binding energy calculated using MMGB/SA was decomposed into energetic contributions from each residue of the binding site ( Figure 7 and Table S2). Binding energy decomposition suggests that all residues of the binding site except for D1532 significantly contributes to the binding affinity of PIIIA to NaV. Binding mode analysis suggests that the D1532 locates in the binding site forming van der Waals contacts with the backbone of the PIIIA (Figure 4). The high desolvation energy of D1532 might result in its unfavorable binding to most of the NaV subtypes, except for NaV1.3 in which the D1532 forms an unstable hydrogen bond with PIIIA O18. An aromatic residue such as Y or F at 401 significantly contributes to the binding affinity of PIIIA to NaV1.1, NaV1.2, NaV1.3, NaV1.4, NaV1.6 and NaV1.7. In our model, the side chain of the aromatic residue at 401 stacks with the side chain of K17, forming favorable van der Waals contacts (Figure 4). By contrast, the C and S merely form few contacts with the side chain of K17, explaining their little energetic contribution to the binding affinity of PIIIA to NaV1.5, NaV1.8 and NaV1.9. A neutral N at 404 is more energetically favorable for the binding of PIIIA to NaV1.x (x = 1, 2, 3, 4, 6, 7) than the positively charged R or K of NaV1.x (x = 5, 8,9). Binding mode analysis indicates that side chain of N can form a hydrogen bond with R14, whereas this hydrogen bond is absent when a positively charged residue is present at 404 (Figure 4). The L and Q at 408 showed only a little difference to the binding affinity of PIIIA to varied subtypes of NaV, despite of their side chain hydrophobicity difference. The L and Q at 408 are at the periphery of the binding site where they both form few van der Waals contacts with the PIIIA explaining their little difference to the binding affinity of PIIIA. The conserved E755 is located at the center of the binding site, substantially contributing to the binding affinity of PIIIA by forming salt bridges with the side chain of R14 with energetic contribution varied from −4 kcal/mol to −7 kcal/mol ( Figure 5). By contrast, the conserved E758 contributes less than E755 to the binding affinity of PIIIA to subtypes of NaV by forming an electrostatic interaction with R12 of PIIIA. A conserved aromatic residue, W761, contributes −2~−1 kcal/mol binding energy to PIIIA by forming van der Waal contacts with side chain of R12. A D at 1241 is conserved among varied subtypes of NaV1 except for NaV1.7, on which an I is present this position. The D1241 locates at the periphery of the binding site where it only contributes −2~0 kcal/mol binding energy to the total binding affinity of PIIIA. The side chain of K17 on PIIIA tilted to the D1241 and could form favorable electrostatic interaction in MD. The D at 1532 is There is no binding or IC 50 data for PIIIA to Na V 1.9. Binding energies calculation from both MMGB/SA and umbrella sampling together predict that PIIIA possesses poor binding affinity to the Na V 1.9.

Binding Energy Decomposition Using MMGB/SA
The binding energy calculated using MMGB/SA was decomposed into energetic contributions from each residue of the binding site ( Figure 7 and Table S2). Binding energy decomposition suggests that all residues of the binding site except for D1532 significantly contributes to the binding affinity of PIIIA to Na V . Binding mode analysis suggests that the D1532 locates in the binding site forming van der Waals contacts with the backbone of the PIIIA (Figure 4). The high desolvation energy of D1532 might result in its unfavorable binding to most of the Na V subtypes, except for Na V 1.3 in which the D1532 forms an unstable hydrogen bond with PIIIA O18. An aromatic residue such as Y or F at 401 significantly contributes to the binding affinity of PIIIA to Na V 1.1, Na V 1.2, Na V 1.3, Na V 1.4, Na V 1.6 and Na V 1.7. In our model, the side chain of the aromatic residue at 401 stacks with the side chain of K17, forming favorable van der Waals contacts (Figure 4). By contrast, the C and S merely form few contacts with the side chain of K17, explaining their little energetic contribution to the binding affinity of PIIIA to Na V 1.5, Na V 1.8 and Na V 1.9. A neutral N at 404 is more energetically favorable for the binding of PIIIA to Na V 1.x (x = 1, 2, 3, 4, 6, 7) than the positively charged R or K of Na V 1.x (x = 5, 8,9). Binding mode analysis indicates that side chain of N can form a hydrogen bond with R14, whereas this hydrogen bond is absent when a positively charged residue is present at 404 (Figure 4). The L and Q at 408 showed only a little difference to the binding affinity of PIIIA to varied subtypes of Na V , despite of their side chain hydrophobicity difference. The L and Q at 408 are at the periphery of the binding site where they both form few van der Waals contacts with the PIIIA explaining their little difference to the binding affinity of PIIIA. The conserved E755 is located at the center of the binding site, substantially contributing to the binding affinity of PIIIA by forming salt bridges with the side chain of R14 with energetic contribution varied from −4 kcal/mol to −7 kcal/mol ( Figure 5). By contrast, the conserved E758 contributes less than E755 to the binding affinity of PIIIA to subtypes of Na V by forming an electrostatic interaction with R12 of PIIIA. A conserved aromatic residue, W761, contributes −2~−1 kcal/mol binding energy to PIIIA by forming van der Waal contacts with side chain of R12. A D at 1241 is conserved among varied subtypes of Na V 1 except for Na V 1.7, on which an I is present this position. The D1241 locates at the periphery of the binding site where it only contributes −2~0 kcal/mol binding energy to the total binding affinity of PIIIA. The side chain of K17 on PIIIA tilted to the D1241 and could form favorable electrostatic interaction in MD. The D at 1532 is conserved across all Na V subtypes. The D1532 contributes little or even unfavorable binding energy to the binding affinity of PIIIA to different subtypes of Na V 1 except for to the Na V 1.3 in which D1532 can form a hydrogen bond with the hydroxyl group of O18. Residues at 1536 are not conserved across different subtypes of Na V . An N at 1536 of Na V 1.4 can form a hydrogen bond with the hydroxyl group of O18, whereas this hydrogen bond is absent when A, S or L are present at the same position. The side chain of L1536 at Na V 1.6 stacks with the side chain of R20 as well as the backbone of PIIIA, and significantly contributes to its binding affinity. A S at 1536 can form hydrogen bond with R20 and significantly contributes to its binding to Na V 1.8, whereas such pairwise interaction is absent at the Na V 1.9. Thus, the same residue at the same positions of different subtypes of Na V might play different roles to the binding affinity of PIIIA. The D at 1541 is conserved across all Na V subtypes from Na V 1.1 to Na V 1.8, whereas an N is at the same position of Na V 1.9. Either an D or an N at 1541 of Na V subtypes contributes to about −2 kcal/mol energy to the binding affinity of PIIIA by forming salt bridges or hydrogen bonds with R20 of PIIIA. conserved across all NaV subtypes. The D1532 contributes little or even unfavorable binding energy to the binding affinity of PIIIA to different subtypes of NaV1 except for to the NaV1.3 in which D1532 can form a hydrogen bond with the hydroxyl group of O18. Residues at 1536 are not conserved across different subtypes of NaV. An N at 1536 of NaV1.4 can form a hydrogen bond with the hydroxyl group of O18, whereas this hydrogen bond is absent when A, S or L are present at the same position. The side chain of L1536 at NaV1.6 stacks with the side chain of R20 as well as the backbone of PIIIA, and significantly contributes to its binding affinity. A S at 1536 can form hydrogen bond with R20 and significantly contributes to its binding to NaV1.8, whereas such pairwise interaction is absent at the NaV1.9. Thus, the same residue at the same positions of different subtypes of NaV might play different roles to the binding affinity of PIIIA. The D at 1541 is conserved across all NaV subtypes from NaV1.1 to NaV1.8, whereas an N is at the same position of NaV1.9. Either an D or an N at 1541 of NaV subtypes contributes to about −2 kcal/mol energy to the binding affinity of PIIIA by forming salt bridges or hydrogen bonds with R20 of PIIIA. Overall, the binding site of NaV is rich of negatively charged residues such as D and E, and thus the binding site shows an overall negative electrostatic potential ( Figure S4). By contrast, the PIIIA contains six positively charged residues and shows the positive electrostatic potential that is complementary to that of the NaV ( Figure S4). Thus, the long range electrostatic interaction can be the driving force for the binding of PIIIA to the NaV at long range distance. At short range distance, the pairwise electrostatic interactions and hydrogen bond interactions contribute most to the binding affinity of PIIIA to NaV subtypes. Binding energy decomposition clearly shows that the conserved residues and the non-conserved residues have synergistic effects to affect the binding affinity of PIIIA to NaV subtypes. Residue side chain energetic contribution to the binding affinity of PIIIA to varied NaV subtypes from NaV1.1 to NaV1.9 are shown in different colors.

Conclusions
A computational method combining homology modeling, docking, molecular dynamics simulations and binding energy calculations was established and successfully applied to determine the specificity of PIIIA to varied NaV subtypes from NaV1.1 to NaV1.9. The predicted binding free energies using umbrella sampling and MMGB/SA were correlated with experimental values giving a correlation coefficient (R 2 ) of 0.68 and 0.41, respectively. The performance of umbrella sampling is Overall, the binding site of Na V is rich of negatively charged residues such as D and E, and thus the binding site shows an overall negative electrostatic potential ( Figure S4). By contrast, the PIIIA contains six positively charged residues and shows the positive electrostatic potential that is complementary to that of the Na V ( Figure S4). Thus, the long range electrostatic interaction can be the driving force for the binding of PIIIA to the Na V at long range distance. At short range distance, the pairwise electrostatic interactions and hydrogen bond interactions contribute most to the binding affinity of PIIIA to Na V subtypes. Binding energy decomposition clearly shows that the conserved residues and the non-conserved residues have synergistic effects to affect the binding affinity of PIIIA to Na V subtypes.

Conclusions
A computational method combining homology modeling, docking, molecular dynamics simulations and binding energy calculations was established and successfully applied to determine the specificity of PIIIA to varied Na V subtypes from Na V 1.1 to Na V 1.9. The predicted binding free energies using umbrella sampling and MMGB/SA were correlated with experimental values giving a correlation coefficient (R 2 ) of 0.68 and 0.41, respectively. The performance of umbrella sampling is superior to MMGB/SA due to the more powerful conformation sampling of the former. Interestingly, the predicted binding affinities given by umbrella sampling is very similar to the experimental values.
Binding modes analysis identified nonconserved residues existing at the interface between PIIIA and Na V 1.x (x = 1 to 9). These nonconserved residues formed nonconserved interactions at the interface between PIIIA and Na V subtypes, which might confer the specificity of PIIIA. Indeed, results from binding energy decomposition showed that these nonconserved residues contributed varied binding energies to the binding affinity of PIIIA to varied Na V subtypes. Intriguingly, we also found that even the conserved residues could contribute varied binding energies to the binding of PIIIA to varied Na V subtypes. Therefore, the conserved and nonconserved residues might play synergistic effects together to influence the specificity of PIIIA to varied Na V subtypes. Overall, results from this study are valuable for our understanding the molecular determinants that confer the specificity of PIIIA to varied Na V subtypes and for further engineering the specificity of the µ-conotoxins to the interested Na V subtypes.