Elucidating molecular mechanisms of protoxin-II state-specific binding to the human Na_V1.7 channel

Voltage-gated sodium (Na_V) channels are responsible for electric signaling and play a key role in various physiological processes (Ahern et al., 2016; Bennett et al., 2019). Na_V channels are involved in the functioning of the central and peripheral nervous systems and the contraction of skeletal and cardiac muscles (Bennett et al., 2019). There are nine known members in the Na_V family, namely Na_V1.1–Na_V1.9 (Bennett et al., 2019). A pore-forming α subunit of a eukaryotic Na_V channel is composed of four homologous domains (DI–DIV), each containing six transmembrane segments (S1–S6), with S1–S4 serving as the voltage-sensing domain (VSD). The S4 segment contains four or five positively charged amino acid residues known as gating charge residues and that act together as a voltage sensor. The S5 and S6 segments, intervening loops, and pore helix form the channel pore containing the ion selectivity filter (SF) region. Membrane depolarization induces the charged residues in the S4 segment to move outward, a movement that is transmitted to the pore domain through the S4–S5 linker leading to channel opening and a rapid influx of sodium ions into the cell. This sodium influx further depolarizes the cell membrane and leads to the generation of an action potential. Na_V channel opening depends on the sequential activation of VSD I through III (Capes et al., 2013; Clairfeuille et al., 2019). In contrast, the activation of VSD IV is coupled with the rapid inactivation of the channel through the release of the isoleucine–phenylalanine–methionine (IFM) or similar motif in the DIII–DIV intracellular loop (inactivation gate; Capes et al., 2013; Clairfeuille et al., 2019; Chanda et al., 2004; Chanda and Bezanilla, 2002). This inactivation process is crucial to stop the influx of sodium ions and enable the repolarization of the membrane to terminate the action potential.

Na_V channel dysfunction arising from mutations can lead to various pathophysiological conditions, including arrhythmias, epilepsy, chronic pain, and insensitivity to pain (Chowdhury and Chanda, 2019). Na_V1.3, Na_V1.7, Na_V1.8, and Na_V1.9 channels have been implicated in pain signaling (Bennett et al., 2019). Transgenic mice lacking Na_V1.7- and Na_V1.8-positive nociceptors have been shown in experiments to exhibit increased mechanical and thermal pain thresholds (Nassar et al., 2004). In particular, the mice showed reduced or no response to inflammatory pain responses evoked by various stimuli (Nassar et al., 2004). Moreover, humans born with loss-of-function mutations in the SCN9A gene encoding for the Na_V1.7 α subunit were discovered to have congenital insensitivity to pain (Drenth and Waxman, 2007). In contrast, gain-of-function mutations have been associated with aberrant chronic pain sensation and related disorders such as inherited erythromelalgia and paroxysmal extreme pain disorder (Golshani et al., 2014; Dib-Hajj et al., 2013). These observations further emphasize the importance of the Na_V1.7 channels to pain perception.

Local anesthetics, such as lidocaine, reduce the transmission of pain signals to the central nervous system by non-selectively blocking various Na_V subtypes (Yang et al., 2020). Due to their lack of Na_V selectivity, side effects are often reported concerning the loss of other sensations, such as the sense of touch. The use of opioid analgesics is associated with side effects such as respiratory depression, constipation, and the development of dependence (Hsu, 2017). Therefore, developing drugs that selectively inhibit Na_V1.7 function may result in strong analgesia without undesirable side effects (Nguyen et al., 2022; Nguyen and Yarov-Yarovoy, 2022). Na_V channels implicated in pain signaling have been discovered to be targets of natural toxins found in species such as spiders and scorpions (Chowdhury and Chanda, 2019). The toxins modulate channel activity by blocking the ion permeation pathway (i.e., pore blockers) or binding to the VSD to alter channel gating kinetics (i.e., gating modifier toxins; Tran et al., 2020).

Protoxin-II (PTx2), a 30-residue peptide (see Data S1) derived from the venom of the Peruvian Green Velvet tarantula (Thrixopelma pruriens), is a gating modifier toxin and has a moderate selectivity for Na_V1.7 versus other Na_V subtypes, making it a suitable scaffold for a peptide-based pain therapeutic (Shen et al., 2019). This toxin interferes with the activation of Na_V channels by binding to the S3–S4 loop in the deactivated VSD II and shifting the voltage dependence of activation to a more positive potential (Shen et al., 2019). PTx2 is 100-fold more potent for Na_V1.7 (IC₅₀ = 0.3–1.0 nM) versus other Na_V subtypes (Nguyen et al., 2022; Schmalhofer et al., 2008; Xiao et al., 2010; de Lera Ruiz and Kraus, 2015). Additionally, PTx2 binds to the deactivated hNa_V1.7 VSD IV inhibiting fast inactivation and producing sustained sodium currents, similar to the effect of α-scorpion toxins (Clairfeuille et al., 2019), but with lower potency (IC₅₀ = 0.24 μM; Xiao et al., 2010). However, this effect of PTx2 is typically masked by the preferential inhibition of Na_V channel activation due to its higher affinity binding to VSD II (Xiao et al., 2010).

By using peptide toxins as a scaffold for designing novel ion channel modulators, researchers can potentially develop more effective and targeted therapies for chronic pain (Nguyen and Yarov-Yarovoy, 2022; Wulff et al., 2019). Potential strategies for pain therapy via hNa_V1.7 gating inhibition include the trapping of VSD II in the deactivated state (Xu et al., 2019; to prevent the channel from activating) or the trapping of VSD IV in the activated state (Ahuja et al., 2015; to keep the channel in the inactivated, non-conducting state). However, current structural studies (Shen et al., 2019; Clairfeuille et al., 2019) have yet to result in a complete atomic-resolution understanding of how PTx2 interacts with Na_V1.7 VSDs in a state-specific manner and with the surrounding lipids. Recently, experimental structures of PTx2 bound to Na_V1.7 VSD II in various states have been resolved using chimeric constructs, grafting parts of the human Na_V1.7 channel sequence into a bacterial Na_V channel (Xu et al., 2019). These studies resulted in X-ray and cryogenic electron-microscopy (cryo-EM) structures of the PTx2-hNa_V1.7-VSD II-Na_VAb complexes in activated and deactivated states (Xu et al., 2019). Yet, these structures do not fully represent a human Na_V channel. Notably, Shen et al. (2019) obtained the first cryo-EM structure of the full hNa_V1.7, resolving structures of PTx2 bound to activated VSD II and VSD IV. However, due to the relatively low resolution of PTx2 densities, these structures did not allow complete atomic reconstruction of the peptide toxin. In addition, the studies did not provide energetic insights into toxin channel and lipid interactions, which are crucial to designing improved peptides as ion-channel modulators.

Computational methods can provide structural insights into the atomic-level mechanisms involved in ion channel function, modulation (Flood et al., 2019), and toxin binding (Henriques et al, 2016; Deplazes, 2017; Mateos and Yarov-Yarovoy, 2023). Moreover, they are useful tools to guide the rational design process toward enhanced peptide variants (Nguyen et al., 2022; Katz et al., 2021). Our study demonstrated that computational structural modeling and molecular dynamics simulations could accurately capture the atomic-resolution molecular mechanisms by which PTx2 binds to hNa_V1.7 VSD II and IV in the activated and deactivated states. Our molecular modeling and simulation insights will be useful to inform future peptide design studies and help researchers identify potential peptide mutations to improve PTx2 selectivity and potency for hNa_V1.7. Moreover, our methodology can be expanded to study other toxin–ion channel interactions to design potential peptide-based therapeutics for diseases such as cancer (Angus and Ruben, 2019), cardiac arrhythmia (Bajaj and Han, 2019; Borrego et al., 2021), epilepsy (Chow et al., 2019; Zhang et al., 2019), and neuropathic pain (Bajaj and Han, 2019; Nguyen et al., 2022; Wulff et al., 2019).

In this study, we used computational modeling methods guided by experimental structural data to investigate molecular mechanisms of PTx2 interaction with the hNa_V1.7 VSD II and VSD IV (see Data S1 for amino acid sequences) in multiple conformational states (Fig. 1). Our primary objective was to generate models of PTx2 bound to VSD II and VSD IV in deactivated and activated conformations as input for molecular dynamics (MD) simulations.

To achieve this goal, we first used cryo-EM structures of PTx2 bound to hNa_V1.7 to both the activated VSD II and VSD IV (PDB accession no. 6J8J; Shen et al., 2019). However, the resolution of the toxin EM density, at ∼5 Å, did not permit complete atomic structure reconstruction in the experimental structure (Shen et al., 2019). To overcome this limitation, we utilized the Rosetta structural modeling software suite (Rohl et al., 2004) and the PTx2–hNa_V1.7 complex electron-microscopy (EM) density data to model the interactions of PTx2 with activated hNa_V1.7 VSD II and VSD IV. For VSD II, we used Rosetta Loop Modeling (Mandell et al., 2009; Stein and Kortemme, 2013) for modeling missing residues in the S3–S4 loop which form receptor site for binding of PTx2 based on experimental structural data for both VSD II (Shen et al., 2019) and PTx2 (Xu et al., 2019).

The PTx2 docking protocol comprised a first step using RosettaDock (Gray et al., 2003; Wang et al., 2007) with the Rosetta membrane energy function (Alford et al., 2020). RosettaDock is a Monte Carlo-based multiscale algorithm that optimizes both rigid-body orientation and side-chain conformations of the docked protein partners. In this step, the ion channel VSD’s structural model was transformed into membrane coordinates by superimposing it into a reference hNa_V1.7 experimental structure downloaded from the Protein Data Bank of Transmembrane Proteins (PDBTM) database (Kozma et al., 2013), which stores experimental structures of membrane proteins in membrane coordinates. The toxin was placed manually in three different initial locations around the binding site in the corresponding VSD using the toxin’s EM density as a reference. The structure of PTx2 was obtained from the PTx2-hNa_V1.7-VSD II-Na_VAb X-ray structure (PDB accession no. 6N4I) (Xu et al., 2019). Subsequently, 20,000 docked models per input were generated with RosettaDock making a total of 60,000 models. The top 10% of models based on Rosetta total_score (total computed Rosetta energy of the protein) were extracted, and the top 100 interface_score (Rosetta interaction energy) models were used as input for the next step. The second step involved EM density fitting of top-docked models to the experimental EM map using Rosetta EM density refinement (Wang et al., 2016b). During EM density refinement, toxin backbone and sidechain conformations are optimized using a modified energy function that accounts for the agreement between the model and the provided experimental EM map. In this step, 500 models were generated per input, making a total of 50,000 models. The top 10% of models based on total_score was initially extracted. The final models were selected based on the recalculated interface_score after EM fitting and on the elec_dens_fast score term (agreement between the Rosetta model and the experimental EM maps).

Furthermore, we docked PTx2 onto the structure of hNa_V1.7 VSD IV trapped in a deactivated state by an α-scorpion toxin (PDB accession no. 6NT4) (Clairfeuille et al., 2019), thereby constructing a model of PTx2 bound to deactivated hNa_V1.7 VSD IV. Note that this experimental structure comprises a human–cockroach chimeric Na_V channel, in which the VSD IV has a human sequence while the remaining part is derived from the cockroach Na_VPaS channel. To model the interactions between PTx2 and VSD IV, the Na_VPaS portion was removed, leaving only the VSD IV. The docking protocol was similar to the two cases above, but without any EM density-related refinement or scoring, since no experimental map is available for this interaction. Finally, we employed the cryo-EM structure of the chimera channel hNa_V1.7-VSD II-Na_VAb (PDB accession no. 6N4R) (Xu et al., 2019), which has the VSD II in a deactivated state, as a template for RosettaCM (Song et al., 2013) homology modeling. This allowed us to generate a model of the deactivated VSD II with the entire hNa_V1.7 VSD II sequence (sequence identity between template and target: 63.7%). We generated 5,000 RosettaCM models, extracted the top 10% based on total_score, and selected the final model based on both low root-mean-square deviation (RMSD) for the template and low total_score. Then, we used RosettaDock (Gray et al., 2003; Wang et al., 2007) with membrane energy function (Alford et al., 2020) to dock PTx2 to the deactivated VSD II of the hNa_V1.7 model described above.

Using the CHARMM-GUI web server (Jo et al., 2008), the VSD-bound-toxin complexes were inserted into tetragonal phospholipid bilayer patches composed of 160 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) molecules in total. The systems were then solvated by a 0.15 M aqueous NaCl solution to mimic the physiological extracellular environment, resulting in molecular systems of ∼63,000 atoms as illustrated in Fig. 2. The protonation state of each residue was assigned under neutral pH conditions. Standard N- and C-termini were set for PTx2 and VSD II/IV. For PTx2, disulfide bonds were added between cysteine residues C2 and C16, C9 and C21, and lastly, C15 and C25.

All-atom MD simulations were performed using the Amber20 software package (Salomon-Ferrer et al., 2013) and the standard all-atom Chemistry at Harvard Macromolecular Mechanics (CHARMM) force fields for proteins (CHARMM36m), lipids (C36), and ions (Huang and MacKerell, 2013; Klauda et al., 2010; Best et al., 2012) as well as the TIP3P water model (Jorgensen et al., 1983). During system equilibration, a 1 kcal/mol/Ų harmonic restraint was applied to all atoms in the system and was gradually reduced to 0.1 kcal/mol/Ų over a period of 90 ns. Subsequently, a production run of 910 ns was conducted, resulting in a total simulation time of 1 μs for each system. All MD simulations were run in the isobaric–isothermal or NPT ensemble at 310.15 K and 1 atm pressure, which was maintained using the Langevin temperature equilibration scheme and Berendsen barostat. Non-bonded interactions were computed up to a 9 Å cutoff. Long-range electrostatic interactions were computed using the particle mesh Ewald (PME) method (Darden et al., 1993), while no long-range correction was applied to van der Waals interactions as suggested for the C36 lipid force field. All covalent bonds involving hydrogen atoms were constrained using the SHAKE algorithm (Ryckaert et al., 1977), allowing the use of a 2-fs time step. The simulation protocol was replicated three times, each with different initial lipid positions and starting velocities of atoms. The replication accounts for potential variability in the lipid arrangement and system dynamics.

Afterward, protein–ligand interactions were characterized in each set of simulations using in-house Python scripts incorporating the protein–ligand interaction profiler (PLIP) Python module (Salentin et al., 2015). Binding free energy calculations were performed using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) approach (Srinivasan et al., 1998; Kollman et al., 2000; Wang et al., 2016a) with all-atom MD simulation trajectories by MMPBSA.py program in Amber Tools (Miller et al., 2012). MM/PBSA was selected due to its computational efficiency, which facilitates the generation of qualitative insights into the energetic components of binding. Moreover, it enables the incorporation of an implicit membrane in its calculations. This feature is particularly advantageous as it fosters a more accurate approximation of the energetics underlying the interaction between PTx2 and hNa_V1.7 embedded in the membrane.

The Chamber module of ParmEd program was used to convert CHARMM-style forcefield files to Amber-style forcefield files without changes in the force field parameters (Shirts et al., 2017). The aqueous solution (with ionic strength of 0.15 M) and lipid membrane were treated implicitly using dielectric constants, denoted ε (water ε_w = 80, lipid bilayer ε_l = 2, and protein ε_p = 4). The solvent probe radius was set to 1.4 Å and the atomic radii were set according to the converted force field parameters.

To obtain the solvation-free energy and the gas-phase electrostatic energy contributions without entropy, the particle–particle particle–mesh (P3M) procedure was used (Lu and Luo, 2003). These calculations were performed with an implicit membrane, where the electrostatic energy includes both reaction field and Coulombic electrostatic energies. Entropy was calculated separately by the interaction entropy method (Duan et al., 2016), where interaction energy between PTx2 and hNa_V1.7 VSD II/IV, including only Coulombic electrostatic and van der Waals energies was needed. To obtain the Coulombic energy separated from the reaction field energy, each system energy was recalculated using the dielectric boundary surface charges method in the implicit ionic solution.

To calculate the energy contribution of each residue to the binding process for each system, we computed the electrostatic and van der Waals interaction energies of PTx2 with POPC lipids and VSD II/IV using AMBER linear interaction energy analysis (Salomon-Ferrer et al., 2013; Roe and Cheatham, 2013)

Molecular graphics visualization was performed using UCSF ChimeraX (Pettersen et al., 2021) to depict all resulting models and illustrate toxin–VSD protein residue interactions. The structural comparison was carried out through the superposition of the structures using the MatchMaker tool within ChimeraX.

Figs. S1, S2, S3, S4, S6, S7, and S8 detail specific PTx2–hNa_V1.7 VSD II, IV, and lipid inter-residue interactions including their frequencies and interaction energies in different states. Figs. S9, S11, and S12 detail a general overview of PTx2–hNa_V1.7 activated VSD II, IV, and lipids inter-residue interactions, now with the VSDs as a part of the full hNa_V1.7 channel. Fig. S10 provides RMSD time series of PTx2 backbone atoms relative to the VSD compared with the initial structures for all the MD simulation systems. Fig. S5 demonstrates time series of MM/PBSA enthalpies and interaction energies for different PTx2–hNa_V1.7 isolated VSD II and IV systems. Fig. S13 provides average MM/PBSA interaction enthalpies and free energies as well as time series of MM/PBSA enthalpies and interaction energies for different PTx2-full channel hNa_V1.7 systems with activated VSD II and IV. Videos 1 and 2 demonstrates structural transitions of PTx2-hNa_V1.7 VSD II and VSD IV complexes between deactivated and activated VSD states using RosettaDock models, respectively. Data S1 includes PTx2, hNa_V1.7 VSD II, and VSD IV amino acid sequences used in this study. Table S1 provides the list of criteria used for detecting different types of non-bonded interactions.

In this study, we employed computational modeling, informed by experimental structural data, to investigate the molecular interactions between PTx2 and hNa_V1.7’s VSDs II and IV in different conformational states (Fig. 1). It should be noted that our focus was solely on individual VSDs and not the whole channel α subunit, an approach selected to optimize the use of computational resources and consistently use available structural information. While experimental structures of activated/deactivated VSD IV and activated VSD II of human Na_V1.7 are available (Shen et al., 2019; Xu et al., 2019), the experimental structure of deactivated VSD II of hNa_V1.7 has yet to be determined. To acquire the necessary structural information, we sought to obtain structures of both hNa_V1.7 VSD II and VSD IV in both activated and deactivated conformations. To obtain a structure of deactivated VSD II, we used the experimental structure of the hNa_V1.7-VSD II-Na_VAb (Xu et al., 2019) as a template for homology modeling (sequence identity: 63.7%) using the RosettaCM (Song et al., 2013). We then compared the positions of gating charges in each VSD by calculating the distances between equivalent Cα atoms after VSD superimposition (Fig. 1). Our analysis revealed that gating charges of VSD II move ∼8.5 Å, on average, between the deactivated and activated states. In the deactivated state, only the first gating charge (R1) is above the hydrophobic constriction site (HCS), while in the activated state R1, R2, and R3 are above the HCS. VSD IV showed a more dramatic displacement of the gating charges (∼13.8 Å) between the activated and deactivated states. Just as in VSD II, in the deactivated conformation, only R1 is located above the HCS; however, in the activated conformation of VSD IV R4 also reaches above the HCS, which contrasts with the activated VSD II structure where only R1 to R3 lies above the HCS. These results suggest that the VSD II in the Na_V1.7 cryo-EM structure (PDB accession no. 6J8J; Shen et al., 2019) is not fully activated, which might be explained by the presence of PTx2 or the lack of membrane voltage during the structure determination process.

Subsequently, structural models of PTx2 in complex with VSD II and VSD IV were generated in both deactivated and activated states as starting points in MD simulations. The methodology utilized for generating these models is outlined in Fig. 3 and briefly summarized as follows: RosettaDock with a membrane energy function was employed for docking PTx2 to the VSDs. The activated conformations of the VSDs (PDB accession no. 6J8J; Shen et al., 2019) with bound PTx2 were further refined using the available experimental EM density maps.

Previously, it was determined that PTx2 interferes with Na_V1.7 activation by binding potently (IC₅₀ = 0.3–1.0 nM) to the S3–S4 loop in VSD II and shifting the voltage dependence of activation to more depolarized potentials (Nguyen et al., 2022; Shen et al., 2019; Schmalhofer et al., 2008; Xiao et al., 2010; de Lera Ruiz and Kraus, 2015). Analysis of the structural models generated here provided important insights into the conformational-dependent binding of PTx2 to Na_V1.7 VSD II and VSD IV. PTx2 binds to the VSD II at the “LFLAD” motif, which comprises residues L812 to D816 and is located within the S3 segment and S3–S4 loop (Xu et al., 2019). The toxin-binding interface can be divided into two distinct regions: polar and hydrophobic. In the PTx2-deactivated VSD II model (Fig. 4 A, panel i), the first gating charge (R1) is observed to be located below the channel-toxin polar binding interface. The electronegative pocket created by VSD II residues E811 and D816 forms hydrogen bonds and salt bridges with the positively charged toxin residues K26 and R22 (Fig. 4 A, panel iii). The hydrophobic interface (Fig. 4 A, panel ii) is composed of M6 and tryptophan residues W5, W24, and W30 (not shown) from the toxin and L812 and F813 from the channel’s VSD II S3 segment. Additionally, W24 was observed to be positioned into the VSD’s top wedge, making hydrophobic contacts with channel residues A766 and L770 from VSD II segment S2 (Fig. 4 A, panel iv).

In the PTx2 - activated VSD II model (Fig. 4 B), the upward movement of the S4 segment and subsequent relocation of the R1 gating charge residue disturbs the electronegative pocket, and the toxin experiences a global shift that pulls it away from the VSD (Video 1). The main interaction in our PTx2-activated VSD II model at the polar interface is between the toxin’s R22 and K26 residues and the channel’s D816 residue (Fig. 4 B, panel iii). Also, the channel’s residue E818 interacts with K28 from the toxin (Fig. 4 B, panel i), although this interaction is not present in any of the experimental structures. The hydrophobic interface (Fig. 4 B, panel ii) comprises π-stacking interactions between the toxin’s W5 and W30 residues and the channel’s F813 residue (Fig. 4 B, panel iv). The overall shift of the toxin resulted in the toxin’s residue W24 no longer interacting with the channel’s residues of the VSD II S2 segment.

Experimental results show that PTx2 can also bind to Na_V1.7 VSD IV and modulate channel inactivation (Xiao et al., 2010), although with much lower potency (IC₅₀ = 0.24 μM). Our models of PTx2 bound to both deactivated and activated VSD IV presented binding interfaces, which can also be separated into a polar interface and a hydrophobic interface. Our models revealed that the position and interface of PTx2 changed significantly more between the deactivated and activated states of VSD IV compared with VSD II (Video 2).

Our model of PTx2 bound to the deactivated VSD IV was built using the experimental structure of the hNa_V1.7 VSD IV in a deactivated state induced by binding of an α-scorpion toxin (PDB accession no. 6NT4) (Clairfeuille et al., 2019), and it is important to note that the actual VSD IV deactivated state to which PTx2 binds may differ. In our final top Rosetta model, which was selected for subsequent MD simulations, PTx2 is embedded in the membrane between the VSD IV S3 and S2 segments. hNa_V1.7 VSD IV residues D1586 and E1589 create an electronegative pocket that accommodates R22 and K26 from the toxin (Fig. 5 A, panel i). Toxin’s R22 extends into the VSD and establishes hydrogen bonds with the channel’s residue D1586, while K26 interacts with G1581 at the top of the S3 segment (Fig. 5 A, panel iii). The hydrophobic interface is more complex in this case (Fig. 5 A, panel ii). PTx2 makes contact with residues in the VSD IV S2 segment, which are mainly hydrophobic interactions mediated by the toxin’s residue W24 and channel’s residues Y1539 and W1538. Finally, toxin’s residues W5, M6, and W30, which also serve as membrane anchors, interact with several hydrophobic residues from the channel located in the VSD IV S3 segment: M1582, F1583, and L1584 (Fig. 5 A, panels ii and iv). In our models, a large conformational rearrangement of the VSD IV occurs between the deactivated and the activated state. As a result, the location of PTx2 also undergoes a dramatic change (Video 2).

In the activated VSD IV state, the toxin is situated at the top of the S3 segment and does not establish interactions with the S2 segment. The polar interface observed in the activated VSD IV state (Fig. 5 B, panel i) is primarily mediated by the interaction between the sidechain of toxin’s residue K26 and the backbone carbonyl oxygen of channel’s S3 segment residue E1589 (Fig. 5 B, panel iii), with PTx2’s residue R22 not appearing to play a significant role in this interaction in contrast to the deactivated state. In the hydrophobic interface (Fig. 5 B, panel ii), the toxin uses its residues W5, M6, and W30 to anchor to the membrane and establish hydrophobic contacts with residues I1588, Y1591, F1592 from the channel’s S3 segment (Fig. 5 B, panel iv). Notably, these hydrophobic channel’s VSD IV residues are positioned far from the binding interface in the deactivated state. In contrast, the channel’s residues previously involved in the hydrophobic interface in the deactivated state (M1582, F1583, and L1584) now lie below the interface and do not interact with the toxin (Fig. 5 B, panel ii).

We conducted MD simulations on the docked PTx2–hNa_V1.7 VSD structures to assess their stability in an environment more closely resembling physiological conditions. These simulations, carried out in three independent 1-µs long replicates, were instrumental in quantifying the energetics of the binding process and investigating the impact of lipids on binding. The simulations showed converging behavior during these timeframes, the details of which are discussed further below. To analyze the simulations, we utilized the atomic coordinates from the entire MD simulation trajectories, excluding 90-ns-long equilibration phases. From this data, contact maps were generated, revealing the residues involved in PTx2–hNa_V1.7 VSD interactions, along with the type (e.g., hydrophobic, hydrogen bond, salt bridges, π-stacking, or cation-π) and duration of these interactions. The criteria used to quantify these interactions are indicated in Table S1. The contact maps are displayed in Figs. S1, S2, S3, and S4.

Throughout the simulation, PTx2 consistently remained bound to VSD II by lodging itself in a cleft located between the S1–S2 and S3–S4 loops, as depicted in Fig. 6, A and D. Hydrophobic contacts constitute the majority of interactions between the toxin and deactivated or activated VSD II (57 ± 4% or 55 ± 5%), followed by hydrogen bond formation (30 ± 2% or 29 ± 6%) and salt bridges (13 ± 2% or 16 ± 8%), as shown in Fig. 6, B and E.

The backbone RMSD profiles, as depicted in Fig. 6, C and F, suggest that the toxin underwent only minor structural changes during the 1-μs simulations, approximately between 1 and 2 Å RMSD. In contrast, the VSDs underwent more considerable structural alterations during the simulations, about 3–4 Å RMSD. The overall RMSD for the complex varied between 3 and 5 Å, with deactivated VSD II systems being on the lower end and the activated VSD II systems being on the higher end. These instabilities could be due to their disconnection from the rest of the ion channel in addition to the flexibility of the loops as well as the pre-S1 helix. However, the RMSD of the toxin-binding region, defined as the area within a 6 Å radius (max distance for interaction detection) of where the toxin interacts with the VSD, experienced only minor deviations from the initial structures, with RMSDs between ∼1 and 2.5 Å.

Figs. S1 and S2 present the contact maps, averaged over three replicas, showcasing the interactions between PTx2 and VSD II in the deactivated and activated states, respectively, along with any interactions with lipids. In both cases, PTx2 formed hydrogen bonds, hydrophobic interactions, and salt bridges primarily with residues on the S3–S4 loop of VSD II (including E811, L812, F813, L814, A815, D816), as depicted in Fig. 7.

In most cases, the RMSD of the PTx2 binding site showed stabilization during the latter half of the simulation (≥500 ns). Based on this observation, we classified contact between two residues as dominant if the average contact duration derived from three replicas, plus standard errors of the mean, exceeds 75% (or 50% for protein–lipid interactions) of the time in this latter half of the simulation. Our subsequent analyses primarily centered on characterizing and visualizing these dominant interactions that drove the binding process. That said, we also paid attention to any noteworthy non-dominant contacts that emerged in the study.

During PTx2 binding to the deactivated VSD II, several key interactions were observed. PTx2 K26 formed a hydrogen bond with VSD II E811 carbonyl located on the S3 helix, while the K27 backbone engaged in hydrogen bonding with F813 carbonyl and A815 amide backbone groups on the VSD II S3–S4 loop (Fig. 7 A). PTx2’s V20 participated in a hydrophobic interaction with VSD II residue A815, and the toxin’s residue R22 formed a salt bridge with VSD II residue D816, as depicted in Fig. 7 C. Notably, as shown in Fig. 7 B, PTx2 was observed to anchor itself to the membrane through hydrophobic interactions involving residue W24 and the surrounding lipids. Additionally, its residue K27 formed salt bridges with nearby lipid head groups. In addition to these dominant interactions, Fig. S1 shows that the toxin’s residues L23 and W24 also sporadically engaged in hydrophobic interactions with VSD II S1–S2 loop residues N763, A766, I767, and L770, effectively creating two anchors that positioned the toxin between the VSD loops.

Intriguingly, when PTx2 binds to the activated VSD II, it exhibits a stronger reliance on interactions with the surrounding lipids rather than VSD II residues. Similar to the previous case, the toxin’s residue K27 backbone engaged in hydrogen bonding with F813 carbonyl and A815 amide backbone groups on the S3–S4 loop (Fig. 7 E). As shown in Fig. 7 D, the PTx2 T8 side chain hydroxyl group engaged in hydrogen bonding, and, to a lesser extent, the K4 side chain amino group formed a salt bridge with nearby POPC phosphate (PO₄) groups. Moreover, as shown in Fig. 7 F, toxin’s residue W30 engaged in hydrophobic interactions with lipids and to a lesser extent, K27 created a salt bridge with POPC PO₄ groups near the VSD II S3 location. At the same time, PTx2 residue W7 participated in cation–π interactions with neighboring lipids close to the VSD II S2 segment (Fig. 7 F). Such interactions contributed to the formation of hydrophobic anchors that helped orient the toxin binding to the VSD II. Although to a lesser extent than the previous case, some hydrophobic interactions were observed involving PTx2 residues W5 and V20 with VSD II residues F813 and A815 (Fig. S2). Additionally, the toxin’s R22 sporadically formed salt bridges with VSD II residues E753 (S1–S2 loop) and D816 (S3–S4 loop).

PTx2 exhibits a dynamic range of binding poses when interacting with VSD IV in the deactivated and activated states. Upon binding the deactivated VSD IV during the simulations, PTx2 progressively adjusted its binding pose, positioning itself deeper into the binding pocket nestled between the S1–S2 and S3–S4 loops as displayed in Fig. 8 A. Hydrophobic contacts formed most of the interactions between PTx2 and the deactivated VSD IV, accounting for 45 ± 5% as demonstrated in Fig. 8 B. These are followed by hydrogen bonds and salt bridges, constituting 32 ± 1% and 22 ± 5% of the interactions respectively. Fig. 8 C depicts considerable fluctuations in the backbone RMSD profiles of the whole complex, the VSD, and the toxin’s binding region throughout the simulations, varying between 2 and 3.5 Å. This significant variability primarily resulted from the flexible S3–S4 loop’s movement as PTx2 embedded itself into the crevice formed between the two VSD loops.

Figs. S3 and S4 present the contact maps, averaged from three replicas, showcasing the interactions between PTx2 and VSD IV in the deactivated and activated states, respectively, along with interactions with lipids. Throughout the MD simulations, PTx2 was bound to the deactivated VSD IV via multiple interactions with the residues on the S3–S4 loop of the VSD. As depicted in Fig. 9 A, PTx2’s K27 formed several hydrogen bonds with F1583 carbonyl and A1585 amide backbone groups, both situated at the VSD IV S3 helix’s terminus where the S3–S4 loop commences. Toxin’s residues K26 and K28 established salt bridges with VSD IV residues E1589 and D1586, respectively, located higher on the S3–S4 loop, as indicated in Fig. 9 B. Additionally, PTx2’s R22 contributed to salt bridge formation with D1586, albeit to a lesser extent than K28 (Fig. S3).

Fig. 9 E illustrates how the toxin bound to the deactivated VSD IV embedded itself within the membrane. PTx2 residues K4 and K27 generated salt bridges with lipid headgroups positioned proximal to the toxin’s binding region on the S3–S4 loop. At the same time, PTx2’s W30 initiated hydrophobic interactions with lipids, effectively functioning as a membrane anchor.

In addition to these dominant interactions, PTx2’s K26 was observed to form a hydrogen bond with VSD IV residue G1581, located lower in the S3 helix (Fig. S3). Moreover, a limited number of interactions of PTx2 with the S1–S2 VSD IV residues were recorded. This includes hydrogen bonding between toxin’s and VSD IV residues T8 and E1534, hydrophobic interactions of L23 and W24 with Y1537 and V1541, respectively, and the formation of a salt bridge between R13 and E1534 as shown in Fig. S3.

Upon binding with the activated VSD IV, as illustrated in Fig. 8 D, the toxin ascended above the membrane, thereby making only a limited number of interactions with lipids. As represented in Fig. 8 E, the interactions established by the toxin with the activated VSD IV were predominantly hydrophobic (49 ± 3%) and hydrogen bonds (46 ± 2%), with salt bridges making up a minor fraction (5 ± 1%). Throughout the MD simulations, the activated VSD IV underwent subtle conformational adjustments, as evidenced by a ∼2–3 Å RMSD relative to the initial structure, shown in Fig. 8 F. However, the toxin’s binding region remained largely unchanged, exhibiting a minor ∼1.5 Å RMSD shift—a fluctuation nearly identical to that of the toxin itself. Interestingly, the RMSD of the complex displayed a greater shift, ranging from 3 to 4.5 Å. Although the structures of the PTx2 and VSD themselves did not exhibit any significant changes, a rotation in the toxin binding position (Fig. 8 D) resulted in the higher complex’s RMSDs observed (Fig. 8 F).

With PTx2 situated atop the activated VSD IV S3–S4 loop, as depicted in Fig. 9 C, the toxin’s K26 side chain amino group formed a hydrogen bond with E1589 backbone carbonyl group of the residue E1589 located on the S3 helix of VSD IV. Concurrently, the toxin’s K27 backbone amide and carbonyl groups established hydrogen bonds with Y1591 carbonyl and V1593 amide backbone groups (Fig. 9 C), respectively, both positioned at the start of the VSD IV S3–S4 loop. Furthermore, PTx2’s V20 engaged in a hydrophobic interaction with VSD IV residue V1593, as visualized in Fig. 9 D. Despite PTx2’s significant elevation above the membrane, which resulted in substantial exposure to the solvent, its C-terminal tail reached down into the lipid bilayer. Here, the toxin’s residue W30 engaged in hydrophobic interactions with nearby lipids in proximity to the VSD IV S3–S4 loop (Fig. 9 F), maintaining sufficient toxin-binding stability. Additionally, marginal hydrophobic contacts involving PTx2’s residue W5 and residue Y1591 on the VSD IV S3-S4 loop were detected as shown in Fig. S4.

We calculated the free energy of binding between PTx2 and different states of hNa_V1.7 VSD II and IV using the implicit-solvent MM/PBSA methodology (Miller et al., 2012) based on our all-atom MD simulation trajectories. We analyzed each component of the free energy by block averages, taken during the latter half of each simulation (500–1,000 ns) when binding stability was typically achieved as assessed by looking at the RMSD profiles (see Figs. 6 and 8), in addition to the enthalpy and the interaction energy time series across all three replicas (Fig. S5). These block averages as well as standard errors of the mean were calculated from all three MD simulation replicas of each system, as illustrated in Fig. 10.

The binding of PTx2 to the activated VSD II and deactivated VSD IV induced the most favorable binding enthalpy (ΔH_MM/PBSA). In contrast, PTx2’s binding to the activated VSD IV, marked by fewer intermolecular interactions in MD simulations, was characterized by the least negative enthalpy. When accounting for the entropic contributions (−TΔS), the resulting binding free energy of PTx2, ΔG_bind, was the most favorable for the deactivated VSD II, with a value of −23.1 ± 0.9 kcal/mol. This was followed by the activated VSD IV (−21.1 ± 1.0 kcal/mol), deactivated VSD IV (−20.1 ± 1.6 kcal/mol), and activated VSD II (−18.5 ± 2.5 kcal/mol) in descending order of favorability. Through paired t tests, we established that a statistically significant difference in ΔG_bind exists (P < 0.05) between PTx2 binding with the deactivated VSD II compared with the activated VSD II and the deactivated VSD IV, respectively.

Our MM/PBSA calculations, showing that PTx2 binds more favorably to the deactivated state of VSD II, align with previous electrophysiology experiments (Xu et al., 2019; Xiao et al., 2010). Those studies demonstrated that PTx2’s antagonism of Na_V1.7 occurs with around 60-fold reduced potency when the membrane potential is held at a voltage that favors VSD activation (Xu et al., 2019). Moreover, another research study suggested that the estimated IC₅₀ for PTx2’s inhibition of Na_V1.7 activation (namely, binding to the deactivated VSD II) is roughly 400-fold lower than the observable IC₅₀ for the inhibition of inactivation (i.e., binding to the deactivated VSD IV; Xiao et al., 2010), corresponding to a ∼3.7 kcal/mol free energy difference in good agreement with our MM/PBSA estimate of 3.0 ± 1.8 kcal/mol.

The time series of MM/PBSA computed enthalpy and interaction energies are shown in Fig. S5. Among the four studied systems, the PTx2 - activated VSD IV system displayed the most stability while the remaining systems showed some degrees of variability across multiple replicas. These variances could stem from the binding mode of PTx2 to the VSD in these cases where the toxin wedged into the space between the two VSD loops. This insertion process caused bending of the S3–S4 loop for the deactivated and activated VSD II as well as deactivated VSD IV systems (as shown in Fig. 6, A and D; and Fig. 8 A), where the attachment occurs, leading to a pattern of contact disruptions and formations. In contrast, the PTx2 - activated VSD IV system displayed the most stability, likely because PTx2 rose above the membrane and engaged in fewer interactions with the activated state of VSD IV (as shown in Fig. 8 D). The limited interactions with the toxin helped preserve the VSD’s structural integrity and kept the binding area largely unaltered.

Then we conducted a linear interaction energy analysis on our all-atom MD simulations using Amber20 (Salomon-Ferrer et al., 2013; Roe and Cheatham, 2013) and the same CHARMM force field parameters as used in the simulations. This allowed us to calculate the electrostatic and van der Waals interaction energies between PTx2 and VSD II/IV residues as well as explicit POPC lipids (not possible with the implicit-membrane MM/PBSA approach), and as depicted in Fig. 11 (for a more exhaustive depiction of specific interactions involving additional residues, refer to Figs. S6, S7, and S8).

Fig. 11 A shows the top 15 PTx2 residues that exhibit the most favorable total interaction energies with activated and deactivated VSD II/IV residues and lipids combined (refer to Fig. S6 to see the full list). The interaction energy for each toxin residue with different VSD states displays a variety of patterns. However, the general trend suggests that most residues have stronger interactions (i.e., more negative interaction energy values) with both the VSD II and VSD IV in their deactivated states compared with their activated states. This could imply that the electrostatic environment of the VSDs in their deactivated state is more favorable for binding those residues. From an energetic perspective, PTx2’s residues R22, K26, K27, K28, and W30 are the main contributors to PTx2–VSD II and IV binding through the formation of hydrogen bonds, hydrophobic interactions, and salt bridges. The most apparent difference is observed in the case of PTx2 binding to the activated state of the VSD IV as the binding pose excludes many toxin residues from forming interactions with the VSD.

The PTx2’s residue K4 also exhibits a significantly more favorable interaction energy when binding to the deactivated state of VSD IV compared with other cases. This can be attributed to the specific binding pose of the toxin in this state, which positions K4 favorably for forming salt bridges with surrounding lipids. The interaction energies with lipids depicted in Fig. S1, S2, S3, and S4 revealed the main PTx2–lipid membrane interaction surface, highlighting the important roles played by toxin residues K4, W5, M6, W7, K27, and W30 in anchoring PTx2 to the cell membrane. The tight binding of PTx2 to the outer leaflet of the cell membrane seems to be critical for the inhibition of PTx2 on hNa_V1.7 channels, which was supported by a mutagenesis study (Henriques et al., 2016) that produced W5Y, W7Y, W24Y, and W30Y PTx2 mutants. In that experimental study, the mutations to tyrosine were performed since tyrosine is more polar than tryptophan and is less suited to bind at the water–lipid bilayer interface (White and Wimley, 1998). These W→Y mutations led to 6- to 291-fold reduction in potency in the inhibition of hNa_V1.7 (Henriques et al., 2016), indicating the importance of these residues for PTx2 binding to hNa_V1.7 in line with our MD simulation predictions.

Despite their significant roles in binding to both states of VSD II and the deactivated state of VSD IV, PTx2 residues K4, W7, T8, R13, R22, and L23 did not form substantial interactions when binding to the activated state of VSD IV as shown in Fig. 11 A. This lack of significant interactions could be attributed to the orientation of the bound toxin, which positioned these residues at a considerable distance from potential binding partners. As a result, their contribution to the binding process is diminished or negligible in this state.

Fig. 11 B provides an overview of the interaction energies between notable VSD II residues with PTx2. When comparing the activated and deactivated states, minor variations are observed in terms of the interaction energies of VSD II residues involved in the binding of PTx2. Notably, residues E811 and D816, located on the S3–S4 loop, exhibit slightly more favorable interaction energies with PTx2 when VSD II is in the deactivated state compared to the activated state. Conversely, residues E753 and E759 on the S1–S2 loop, as well as E818 on the S3–S4 loop, show stronger interactions with PTx2 when VSD II is activated compared with when it is in the deactivated state.

Interestingly, Fig. 11 C, which illustrates the interaction energies of notable VSD IV residues with PTx2, demonstrates a distinct difference in their pattern when PTx2 binds to the activated and deactivated states of VSD IV. This disparity in interaction energies can be attributed to the positioning of PTx2 within the membrane. In the deactivated state (Fig. 8 A), PTx2 binds deeper within the VSD IV cleft, leading to a stronger interaction with the VSD residues. Conversely, in the activated state, PTx2 binds more superficially to VSD IV, resulting in fewer interactions and primarily involving VSD residues located at the top of the S3–S4 loop (Fig. 8 D). Additionally, the analysis of PTx2’s interactions with the lipid membrane reveals that when binding to VSD IV in the activated state, PTx2 exhibits significantly fewer interactions with the lipid membrane compared with the deactivated state (Fig. 9 F versus Fig. 9 E). In conclusion, the variations in interaction energies and the involvement of specific VSD residues, as well as the differential interactions with the lipid membrane, contribute to the distinct binding characteristics of PTx2 to VSD IV in different states.

While the toxin binding regions remained largely stable, structural instability of the rest of the VSDs was observed in certain simulations, likely due to their detachment from the full ion channel’s structure. Control MD simulations with the full hNa_V1.7 channel were conducted to evaluate the toxin–VSD complex stability and the subsequent effect on their binding. Given the unavailability of complete hNa_V1.7 or homologous ion channel structures with deactivated VSDs, our simulations were limited to activated VSD II and IV integrated into the full hNa_V1.7 channel.

As shown in Fig. S9, analysis of three 300-ns full-channel simulation replicas showed improved stability in VSDs, with backbone RMSDs relative to the starting structures of around 2–3 Å for activated VSDs II and IV (see Fig. S9, C and F) versus 3–4 Å for isolated VSDs (see Fig. 6 F and Fig. 8 F). The toxin’s binding region RMSDs stabilized to around 2 Å (VSD II) or 1–1.5 Å (VSD IV) by the end of the MD simulations, comparable to simulations with isolated VSDs.

Fig. S10 offers an analysis of PTx2’s positional variations relative to the VSD in scenarios where the VSD was either isolated or part of the full channel. This was done by calculating the RMSD of PTx2’s backbone atoms when the trajectory was aligned to the VSD. Overall, PTx2’s binding to VSD II and IV in their activated state was characterized by greater fluctuations. This could be attributed to multiple factors. For instance, when PTx2 interacted with the activated VSD II, it engaged more with surrounding lipids than when binding to the deactivated state; the lipids’ inherent mobility could contribute to PTx2’s positional instability. In the case of PTx2 binding to the activated VSD IV, the toxin was observed to rise well above the membrane surface, resulting in fewer VSD contacts and increased exposure to the aqueous environment, which could explain the noted fluctuations. Despite these movements, the RMSD of PTx2 generally stabilized as the simulations progressed, although there were exceptions such as replica 1 of the deactivated VSD IV and replica 2 of the activated VSD IV—both involving isolated VSD systems. Nonetheless, these instances of fluctuation appeared not to compromise the toxin’s binding as the time series of interaction energies for these replicas showed minimal disruption, as shown in Fig. S5, C and D.

Fig. S10, B and D present a comparison of the RMSD time series data for PTx2 in its interaction with both isolated and non-isolated (attached to the full ion channel) activated VSD II and IV systems. When bound to the activated form of VSD II, PTx2 showed only slight positional deviations, regardless of whether the VSD was part of the full ion channel or isolated. In contrast, PTx2’s binding to the activated VSD IV was more stable when the VSD was part of the full channel than in the isolated configuration. This finding is particularly notable given that the activated VSD IV generally displayed more structural rigidity due to its fewer interactions with the toxin in both scenarios. Nonetheless, in the isolated VSD scenario, PTx2 underwent a subtle rotational movement at the binding site (illustrated in Fig. 8 D), a phenomenon not observed when the VSD was integrated within the full channel (see Fig. S9 D).

Further examination of the MD simulation trajectories revealed that when the activated VSD IV was isolated, the toxin’s high exposure to the solvent and minimal VSD binding led to the toxin moving around and causing the isolated VSD’s S3 and S4 helices to sway along (as shown in Fig. 8 F, with an RMSD for the VSD of ∼3 Å). Conversely, when the activated VSD IV was part of the complete channel, the movement of the VSD’s S4 was restricted (as indicated in Fig. S9 F, with an RMSD for the VSD of about 2 Å), resulting in the toxin maintaining its initial binding orientation without significant deviation.

Figs. S11 and S12 depict binding interactions between the toxin and the full-channel activated VSD II and IV during the latter half of the simulations when the protein structures demonstrated stability. These results largely align with those from the isolated VSD simulations. Specifically, for VSD II, PTx2’s K27 hydrogen bonded with VSD’s F813 and A815; W5, V20, and K26 had hydrophobic interactions with VSD S3–S4 loop residues; and L23, W24 with VSD S1–S2 loop residues. PTx2’s R22 and K28 formed salt bridges with VSD’s D816 and E818, respectively. PTx2’s interactions with lipids included T8 via hydrogen bonding, W30 via hydrophobic interactions, K4 and K27 via salt bridges, and W7 via cation–π interactions. In VSD IV, hydrogen bonds included K26-E1589 and K27-Y1591/V1593 pairs; hydrophobic interactions involved PTx2’s W5, V20, K26, and K27 with VSD’s Y1591, F1592, and V1593; and a salt bridge formed between PTx2’s K26 and VSD’s E1589. W30 also interacted with lipids via hydrophobic interactions.

Minor differences were observed compared with the isolated VSD simulations. For example, in the VSD II case, the toxin’s K26 formed a salt bridge with VSD II’s D816, whereas R13 formed a salt bridge and W30 engaged in cation–π interactions with lipids, forming interactions not seen in the isolated VSD case. In the VSD IV case, in addition to W30, W5 also formed hydrophobic interactions with lipids, unlike in the isolated VSD scenario.

Fig. S13 details the calculation of MM/PBSA binding energies for PTx2 with the activated VSD II and IV attached to the hNa_V1.7 channel. Even though the PTx2 interaction with the activated VSD II showed a more favorable binding enthalpy, ΔH, there was no notable distinction in the overall free energy of binding, ΔG, when comparing the activated VSD II with the activated VSD IV. This aligns with the energy computations for the isolated VSD systems shown in Fig. 10, which indicate similar binding free energies of the toxin for both VSD II and IV in the activated state.

In conclusion, integrating the VSDs into the full channel improved stability, although the toxin-binding interaction outcomes largely mirrored those observed in isolated VSD MD simulations. As we do not have the full hNa_V1.7 structure with VSDs in the deactivated state, our study primarily used isolated VSDs, whose structures we have in both the activated and deactivated states, to allow for a consistent comparison.

Using existing experimental data as a guide, we generated models of PTx2 bound to human Na_V1.7 VSD II and IV in both the activated and deactivated states. These models accurately reproduced the crucial toxin–ion channel’s VSD II interactions, as previously elucidated through mutagenesis experiments and structural studies involving chimeric ion channels (Xu et al., 2019; Shen et al., 2019). Our findings also shed new light on PTx2 interactions with VSD IV, unveiling substantial differences compared to its interactions with VSD II.

MD simulations can provide a dynamic perspective on toxin–receptor interactions, transcending the static results from experimental structures and docking. This is especially important when the binding site is surrounded by lipids, making the fluid nature of lipid membranes play a crucial role. As lipids are continually in motion, their interactions with the toxin can significantly influence and even shift the binding of the toxin to the receptor, a phenomenon not captured in static structures or even typical docking runs but evident in MD simulations. Indeed, we discovered several new insights into PTx2–hNa_V1.7 interactions by performing the MD simulations. For example, initially, PTx2’s residue K27 showed no interaction with VSD II in both states and M6 interacted with hydrophobic VSD II residues. However, after 1 µs of the MD simulation, K27 extensively interacted with VSD II, while M6’s interactions diminished. In the case of PTx2 binding to deactivated VSD IV, MD revealed PTx2’s residue K28 as the dominant contact with VSD’s D1586, overshadowing the initially observed PTx2’s R22 interaction. For the activated VSD IV, interactions involving PTx2’s K27 and V20 emerged only in MD simulations. Additionally, MD runs can detail the roles of individual residues in interaction energies, deepening our grasp of the binding mechanics. This detailed perspective, highlighting the influence of dynamic lipids, is vital for understanding molecular interactions and subsequent therapeutic peptide design.

Based on our atomistic structural modeling and MD simulations, we established that PTx2 primarily binds to VSD II in different conformational states through interactions with residues on the S3–S4 loop of the VSD II. In both states, PTx2’s residue K27 formed hydrogen bonds with the channel’s residues F813 and A815 located on the S3–S4 loop. In the deactivated state of VSD II, a hydrogen bond is established between residue K26 of PTx2 and residue E811 of VSD II. On the contrary, when VSD II is activated, K26 interacts with VSD II residue L812 instead, and this interaction was less frequent and more variable. Similarly, PTx2’s R22 established a more robust and consistent salt bridge with D816 (on the S3–S4 loop) when bound to the deactivated VSD II state compared with the activated state, as it occasionally shifted its interaction to E753 (on the S1–S2 loop). Also, PTx2 exhibited a hydrophobic interaction via residue L29 with the channel’s residue L814 in the deactivated state of VSD II, an interaction not observed in the activated state.

PTx2 interacts with VSD IV primarily through residues on the VSD IV S3–S4 loop, but significant differences were observed in the toxin’s binding to the deactivated and activated states of VSD IV. For example, PTx2’s K26 formed a hydrogen bond with G1581 (located in the S3–S4 loop) in the deactivated state and with E1589 (also in the S3–S4 loop) in the activated state of VSD IV. In the activated state, K26’s interaction with E1589 was more frequent and consistent, suggesting stronger binding. Another differential interaction can be seen with the toxin’s residue K27. When binding to the activated VSD IV, K27 interacted with three different residues in the S3–S4 loop, Y1591, V1593, and F1592, while in the deactivated state, it interacted mostly with F1583 and A1585 in the S3–S4 loop only.

A significant advancement highlighted in this study is the novel atomistic understanding of PTx2’s interactions with surrounding POPC membrane lipids upon binding to the Na_V1.7 channel. Lipid interactions play a vital role in stabilizing the position and orientation of PTx2 for binding to VSD II and IV in different conformational states. These novel insights from simulations could open a new avenue for enhancing PTx2 binding to the channel by optimizing the toxin’s interactions with nearby membrane lipids (Henriques et al., 2016).

Upon PTx2’s binding to the activated VSD II, PTx2’s residue T8 formed a hydrogen bond while K4 and K27 formed salt bridges with nearby POPC PO₄ groups. W7, W24, and W30 engaged in hydrophobic interactions with the lipids’ fatty acid tails, and W7 formed cation–π interactions with the POPC choline group. When PTx2 binds to the deactivated VSD II, PTx2’s W5 and W24 engaged in hydrophobic interactions while K4 and K27 formed salt bridges with POPC PO₄ groups. T8 still engaged in hydrogen bonding with nearby lipids, but no hydrophobic interactions originating from W7 and W30 were observed, unlike the previous case. However, PTx2’s R13 was observed to interact sporadically with lipids through salt bridge formation.

When PTx2 binds to the activated VSD IV, the only interaction with lipids was a hydrophobic one originating from W30, which was stronger and more consistent than when PTx2 binds to the deactivated state. However, when PTx2 binds to the deactivated VSD IV, salt bridges with POPC head groups were formed by toxin’s residues K4 and K27 (and with R13 and K28 to a much lesser extent). It is noteworthy that K27–POPC salt bridge interactions were also substantially more frequent and consistent in the deactivated VSD IV than in the deactivated VSD II. Hydrophobic interactions with POPC were observed not only from W30, like the activated VSD IV case, but also from W5 and W24, which are unique to the deactivated VSD IV case.

In summary, this first view into the lipid modulation of toxin binding may constitute a novel and emerging area of research for optimizing peptide toxin’s specificity and selectivity.

Lastly, we complemented the data with specific binding energetic contributions from the individual toxin’s and ion channel’s residues in different VSD conformational states, enabling the identification of potential mutation candidates to enhance or weaken the toxin’s state-specific binding. This information can support the design of toxin analogs with enhanced potency, holding promise as a way to identify candidates for the development of innovative treatments for chronic pain. Our analysis of PTx2 binding energetics to different states of VSD II and IV indicates that PTx2 binds most favorably to the deactivated VSD II, followed by the activated VSD IV. A substantial difference in the free energy of PTx2 binding is observed when comparing the deactivated VSD II to both the activated VSD II and the deactivated VSD IV. PTx2 binding to and trapping VSD II in its deactivated state can prevent the opening of the hNa_V1.7 channel. Likewise, trapping VSD IV in its activated state can keep the channel in a non-conducting, inactivated state. The net effect is a decrease in the Na⁺ current flowing through the Na_V1.7 channels when PTx2 is introduced, which aligns with experimental results (Xiao et al., 2010). From a therapeutic standpoint, optimizing PTx2’s interactions with these states could enhance its effectiveness in pain therapy through hNa_V1.7 channel inhibition. Furthermore, experimental findings demonstrate that the efficacy of PTx2 as an inhibitor of hNa_V1.7 is significantly influenced by its binding to the lipid membrane (Henriques et al., 2016). Therefore, improving PTx2 interactions with cell membrane lipids could also augment PTx2’s inhibition of hNa_V1.7.

Based on our analysis in Fig. 11 A and Fig. S6, we have identified certain residues of PTx2 that interact with either the activated or deactivated states of VSDs and lipids but exhibit relatively weak interaction energies. For residues that primarily interact with membrane lipids (T8, R13, E17), introducing mutations that enhance their hydrogen bonding interactions with lipids for T8, remove the negative charge for E17, or diminish preferential state-specific salt bridge formation for the deactivated VSD IV system could be beneficial. V20 interacts exclusively with the VSDs through hydrophobic interactions but has remarkably low interaction energy despite its relatively long contact duration. Therefore, a non-conservative mutation could be considered. A recent mutagenesis study achieved a more potent and selective peptide by performing the V20R mutation (Nguyen et al., 2022).

A promising strategy entails fine-tuning state-specific interactions. Mutations could be made to increase PTx2’s affinity for states that decrease Na⁺ conduction by preventing channel activation or trapping it in an inactivated state (deactivated VSD II and activated VSD IV, respectively) compared with those that facilitate Na⁺ conduction (activated VSD II and deactivated VSD IV). Residue L23 demonstrated hydrophobic interactions with both the VSDs and the lipid molecules in the membrane. However, it favors binding to the activated VSD II and deactivated VSD IV which are both undesirable. In addition, K4, R13, K26, and K28 exhibit significantly higher binding energies upon binding to the deactivated VSD IV compared with the other states. This heightened affinity is undesirable as it impedes channel inactivation. Consequently, careful mutation of these residues could serve to increase PTx2’s affinity for the desired states of the channel. As an example, a mutagenesis study discovered that K26R and K28E resulted in improved PTx2 potency (Nguyen et al., 2022). Another experimental study also revealed that charge-reversing and charge-neutralizing mutations of PTx2’s basic residues R22 and K26 lead to significant loss of potency (Xu et al., 2019), which was also corroborated by a computer simulation study of PTx2 mutant–hNa_V1.7-Na_VAb-deactivated VSD II system using free energy perturbation (FEP) (Katz et al., 2021). However, this MD simulation method, more rigorous and computationally expensive than a posteriori interaction energy analysis used in this work, was unable to predict gain in potency for some charge-preserving K26 and R22 mutations (e.g., K26R and R22noR, where norR is norarginine) (Katz et al., 2021). This suggests that using PTx2–hNa_V1.7 VSD II and IV systems in different conformational states as we explored in this work might be crucial for a more accurate estimate of peptide toxin mutation effects using FEP or other approaches as will be explored in our follow-up studies.

By strategically introducing mutations to enhance the interactions of these residues, there is a possibility of enhancing the affinity of PTx2 towards the VSDs, their different conformational states, and/or lipids. These mutations can be designed to optimize specific interactions, such as salt bridges, hydrogen bonding, electrostatic, hydrophobic, and/or other types of interactions, involved in PTx2’s binding to the VSDs and lipids as indicated in the contact maps in Figs. S1, S2, S3, and S4. Nonetheless, it is vital to balance this enhancement with the potential effects on state specificity and ion channel subtype selectivity for toxin–receptor binding, as well as the net effect on the protein’s overall structure and stability.

Using computational modeling and simulation tools, we can model the atomic-resolution structures and relative energetics of peptide toxin interactions with hNa_V1.7 in different conformational states. This allows for the accurate prediction of key binding interactions, their subsequent experimental validation, and design of more potent and selective peptides targeting hNa_V1.7 or other ion channels as we did recently (Nguyen et al., 2022).

Although the computational workflow in our study can provide valuable insights into the complexities of peptide–receptor interactions at the atomic resolution, several drawbacks are associated with each methodology. Given the absence of complete hNa_V1.7 channel structures with VSDs in the deactivated state, our modeling and simulations were limited to the isolated VSDs of the channel. This limitation could potentially introduce instabilities during the MD simulations. Nonetheless, we observed only minor displacements of the PTx2’s binding regions during the majority of our simulations as demonstrated by the RMSD profiles shown in Figs. 6 and 8. We also performed preliminary MD simulations of a full hNa_V1.7 channel with VSDs II and IV in the activated state as shown in Figs. S9, S10, S11, S12, and S13, demonstrating fairly similar results to the isolated VSD simulations with somewhat higher stabilities (lower RMSDs) and a few additional interactions observed. In the future, we will expand our study by modeling and simulating full hNa_V1.7 channel structural models in different conformational states.

Peptide backbone flexibility, especially at the N- and C-termini, is difficult to capture during the protein–protein docking and, therefore, RosettaDock might not be able to sample the full range of possible toxin conformations. Additionally, Rosetta scoring uses an implicit membrane model that might not account for energetic contributions involving membrane lipids that could be important to stabilize the toxin’s binding pose. Although a lipid membrane can be explicitly included in MD simulations, simulation accuracy often depends on the quality of the force field and parameters used in the simulation, and it can be difficult to determine the best set of parameters for any given system (Guvench and MacKerell, 2008). Fully atomistic MD simulations are limited by their timescale (ns to μs), and hence enhanced sampling MD simulation techniques such as the string method with the swarm of trajectories (Chen et al., 2022), weighted ensemble (WE; Zuckerman and Chong, 2017), Gaussian accelerated MD (GaMD; Wang and Miao, 2020), or their combination (Ahn et al., 2021) need to be employed to observe peptide–protein binding and large-scale protein conformational transitions.

For estimating the free energy of binding, MM/PBSA methodology may suffer from accuracy problems, especially for predicting absolute binding free energies as was demonstrated by a study done by Sheng et al. (2021). The accuracy can be improved such as by damping the solvated and Coulombic terms for highly charged systems (Sheng et al., 2021; Spiliotopoulos et al., 2016), using residue type-specific dielectric constants (Liu et al., 2019), and applying potentially more accurate entropy calculation methods (Sun et al., 2018). Yet, the MM/PBSA methodology employed here was successfully utilized in our recent study (Han et al., 2023) to correctly predict the conformation-dependent binding of a small molecule ligand norepinephrine and large stimulatory G protein to β-2 adrenergic receptor, another integral membrane protein. To predict the effect of the toxin’s mutations on state-specific ion channel binding, a potentially more accurate but computationally expensive MD simulation technique, FEP (Kollman, 1993), can be used, which was recently employed to study relative energetics of PTx2 mutants–hNa_V1.7-Na_VAb deactivated VSD II interactions (Katz et al., 2021). A similar approach can be used in our future studies focused on optimizing specific toxin–ion channel interactions, while in this work, we aimed to predict molecular determinants of those state-specific interactions using a consistent set of wild-type toxin–hNa_V1.7 VSD II and IV structural models and microsecond-long MD simulations. Other molecular simulation and analysis methodology advancements such as using enhanced sampling techniques GaMD, WE, or their combination (Zuckerman and Chong, 2017; Ahn et al., 2020, 2021; Wang and Miao, 2020) will be explored in our follow-up studies to thoroughly explore both thermodynamics and kinetics of toxin–ion channel interactions with the current work providing a necessary framework and achieving a good agreement with several experimental findings (Xiao et al., 2010; Xu et al., 2019).

In conclusion, our molecular modeling and simulation protocol can be a helpful guide for future designs of more potent therapeutic peptide variants with improved selectivity for different states of the hNa_V1.7 channel and other Na_V channel subtypes as will be explored in follow-up studies. In addition, our results demonstrate that the interactions between PTx2 residues and lipids are important for state-specific binding to hNa_V1.7 VSDs and should be considered in future ion channel inhibitor studies as well. These findings open the possibility of improving the balance of VSD and lipid-binding capability of the toxin and, consequently, its efficacy as an ion channel inhibitor. We anticipate that this methodology can be readily expanded to study the interactions between other peptides and ion channels. Advancing our understanding of the molecular mechanisms of ion channel modulation by peptide toxins will allow us to harness the full potential of these molecules for therapeutic applications.

The data underlying Figs. 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, and 11, as well as Figs. S1, S2, S3, S4, S5, S6, S7, S8, S9, S10, S11, S12, and S13, and Videos 1 and 2 are openly available in the Dryad repository https://doi.org/10.5061/dryad.jq2bvq8gj.

Crina M. Nimigean served as the editor.

We thank all members of the I. Vorobyov, C.E. Clancy, and V. Yarov-Yarovoy laboratories for helpful discussions.

This work was supported by National Institutes of Health (NIH) Common Fund Grant OT2OD026580 (to C.E. Clancy and I. Vorobyov), National Heart, Lung, and Blood Institute (NHLBI) grants R01HL128537, R01HL152681, R01HL085844, and U01HL126273 (to C.E. Clancy, V. Yarov-Yarovoy, and I. Vorobyov), American Heart Association Career Development Award grant 19CDA34770101 (to I. Vorobyov), National Science Foundation travel grant 2032486 (to I. Vorobyov), UC Davis Department of Physiology and Membrane Biology Research Partnership Fund (to C.E. Clancy and I. Vorobyov) as well as UC Davis T32 Predoctoral Training in Basic and Translational Cardiovascular Medicine fellowship supported in part by NHLBI Institutional Training Grant T32HL086350 (to K. Ngo and K.C. Rouen) and UC Davis Chemical Biology Program fellowship supported in part by National Institute of General Medical Sciences Institutional Training Grant 5T32GM136597-02 (to K.C. Rouen). Computer allocations were provided through Extreme Science and Engineering Discovery Environment grant MCB170095 (to I. Vorobyov, C.E. Clancy, and V. Yarov-Yarovoy), Texas Advanced Computing Center Leadership Resource and Pathways Allocations MCB20010 (I. Vorobyov, C.E. Clancy, and V. Yarov-Yarovoy), Oracle for Research fellowship and cloud credits award (to I. Vorobyov, C.E. Clancy), Pittsburgh Supercomputing Center (PSC) Anton 2 allocations PSCA17085P, MCB160089P, PSCA18077P, PSCA17085P, and PSCA16108P (to I. Vorobyov, C.E. Clancy, and V. Yarov-Yarovoy). Anton 2 computer time was provided by the PSC through Grant R01GM116961 from the NIH. The Anton 2 machine at PSC was generously made available by D.E. Shaw Research.

Author contributions: K. Ngo, D. Lopez Mateos, S.-H. Ahn, H. Wulff, C.E. Clancy, V. Yarov-Yarovoy, and I. Vorobyov designed the research; K. Ngo, D. Lopez Mateos, Y. Han, and K.C. Rouen performed the research and analyzed the data. All authors wrote and revised the manuscript and approved its final version.

The data underlying Figs. 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, and 11, as well as Figs. S1, S2, S3, S4, S5, S6, S7, S8, S9, S10, S11, S12, and S13, and Videos 1 and 2 are openly available in the Dryad repository https://doi.org/10.5061/dryad.jq2bvq8gj.