Drugs exhibit diverse binding modes and access routes in the Na_v1.5 cardiac sodium channel pore

Local anesthetic, anti-seizure, and anti-arrhythmic drugs are commonly used to treat a variety of cardiac and nervous system pathologies such as long QT syndrome (Benhorin et al., 2000) and epileptic seizures (Pal et al., 2021). These drugs target voltage-gated sodium (Na_v) channels, which are proteins that form a pore in the membrane to selectively control sodium ion conduction, initiating the action potential in excitable cells. Mammals have nine Na_v subtypes, termed Na_v1.1–Na_v1.9, each with slightly different biophysical characteristics, which are expressed in neurons of the central nervous system, as well as peripherally, in the spinal cord, muscle, and cardiac tissue. Several subtypes, Na_v1.3, Na_v1.5, and Na_v1.7, have now been resolved via cryo-EM with various drug molecules inside the pore (summarized in Table 1). All Na_v channel structures reveal a consistent architecture of peripheral voltage-sensing domains, surrounding a highly conserved central pore. The channel is divided into four homologous domains (DI–DIV), with each domain consisting of six transmembrane helices (S1–S6), the first four of which form separate voltage sensors while the pore is formed from S5 and S6 helices from all four domains (Fig. 1, A, C, and D). The ion-conducting pore contains a narrow selectivity filter (formed by the four re-entrant pore loops), a wider central cavity, and an activation gate formed by the cytoplasmic end of the S6 helices that opens to facilitate sodium ion influx. These structures also reveal four lateral fenestrations, which are orthogonal portals that open from the membrane to the central pore cavity (Fig. 1, B and D).

Small molecule compounds that inhibit Na_v channels can enter and block the pore via two distinct mechanisms: (1) “use-dependent block,” where a drug enters the cavity from the cytoplasm via the open pore gate; or (2) “tonic block,” where a drug diffuses via the membrane into the cavity and does not require gate opening. While first hypothesized 50 years ago (Hille, 1977), experiments (Gamal El-Din et al., 2018) and simulations (Raju et al., 2013; Boiteux et al., 2014; Martin and Corry, 2014) in ancestral bacterial Na_v channels, as well as, the numerous cryo-EM structures, have shown that the fenestrations provide an accessible route for drugs to enter the cavity during tonic block. These drugs block ion conduction by binding in the central cavity, forming interactions with various pore-lining residues, and either physically occluding the pore or stabilizing the non-conductive inactivated state. Notable conserved S6 helix residues are proposed to form the drug receptor site, including F1762 and Y1769 on DIV, and additionally, L1464, N1465, and I1468 on DIII, and I409 on DI (Fig. 1, A and C; numbering based on rat Na_v1.5) (Ragsdale et al., 1994, 1996; Qu et al., 1995; Wang et al., 1995; Li et al., 1999; Yarov-Yarovoy et al., 2001, 2002; Liu et al., 2003). Pore inhibitors typically bind with the greatest affinity to the inactivated state, stabilizing the pore in this conformation and delaying their transition back to the channel’s resting state (Tucker and Mather, 1979; Földi et al., 2021).

Numerous Na_v channel pore blockers are used for different clinical purposes. Lidocaine was one of the first discovered local anesthetics identified to block the pores of Na_v channels. Derivatives of the lidocaine structure, etidocaine, and bupivacaine were shown to have longer-lasting anesthesia due to higher affinity binding and slower dissociation (Tucker and Mather, 1979). Anticonvulsant compounds are typically more hydrophobic than local anesthetics and include early compounds such as phenytoin and carbamazepine, which are frequently used to treat a range of epilepsy types (Nevitt et al., 2019). As compound screening abilities improved, newer anti-convulsant drugs were discovered and approved, such as lacosamide, which had a unique anticonvulsant profile and improved safety (Stöhr et al., 2007). The third type of Na_v inhibitors are antiarrhythmics, which can be classed into 1a, 1b, and 1c depending on their effect on the cardiac action potential (Lei et al., 2018). Class 1a antiarrhythmics, such as quinidine, have intermediate association/dissociation kinetics and prolong the cardiac action potential. Lidocaine is also considered a class 1b anti-arrhythmic, with fast association/dissociation kinetics shortening the cardiac action potential by preventing late sodium currents. Class 1c antiarrhythmics, such as flecainide, have potent antiarrhythmic properties, with high affinity and prolonged duration of action (slow association/dissociation kinetics) compared with lidocaine, but with less effect on individual action potential duration to its long-lasting action. Other small molecule drugs have recently been of interest for their ability to uniquely inhibit Na_v channels, like riluzole, a neuroprotective drug approved for amyotrophic lateral sclerosis (ALS), has recently been shown to have a distinct modulatory effect on Na_v channel function, without fully inhibiting ion conduction (Földi et al., 2021).

Although it is understood that Na_v channel-inhibiting drugs act via accessing and binding in the pore, we still lack a detailed understanding of how these compounds reach the cavity and the precise locations in the pore to which they bind. In particular, it is not clear how such a large variety of compounds can act within the pore, nor the mechanisms by which they yield disparate effects. Recently, several structures of Na_v channels with different inhibitors bound to the pore have been resolved, allowing partial answers to these questions. Structures now exist containing the common inhibitors, flecainide (Jiang et al., 2020) and quinidine (Li et al., 2021), in Na_v1.5, as well as, lamotrigine, riluzole (Huang et al., 2023b), bupivacaine, lacosamide (Wu et al., 2023), and cannabidiol (Huang et al., 2023a) in Na_v1.7, which show a diversity of poses that are occupied within the pore cavity.

Computational studies have also aided in understanding aspects of drug binding in the Na_v channel pore that are hard to capture in structural studies (Chen et al., 2017). Because of the slow timescale of drug binding relative to typical molecular dynamics simulations, these often use enhanced sampling methods to better sample binding positions in the pore. Docking studies of pore blockers showed that compounds have diverse poses but share some common interactions, such as with an aromatic phenylalanine, or with/replacing sodium ions in the selectivity filter (Tikhonov and Zhorov, 2017). Simulations on the early bacterial structures examined the fenestration size (Kaczmarski and Corry, 2014) and the ability of compounds to pass through the fenestrations and potential binding locations using the approach of umbrella sampling (Martin and Corry, 2014) and flooding with high concentrations of drug (Raju et al., 2013; Barber et al., 2014; Boiteux et al., 2014; Nguyen et al., 2019) and metadynamics (Smith and Corry, 2016). Using replica exchange simulations, we have previously shown that different inhibitors share a set of common binding sites within the pore cavity but each compound is seen to fluctuate in occupying several diverse sites (Buyan et al., 2018, 2021). We were also able to show differences between neutral and charged inhibitors, with the latter interacting more strongly with the oppositely charged selectivity filter region. More recently, we used metadynamics to quantify the energetics of lidocaine inside the pore (Tao and Corry, 2022). We showed that lidocaine can enter through specific fenestrations and occupy multiple equally favorable binding poses within the pore. This further supports the conclusion that many drugs can enter the pore through the fenestrations and demonstrates the utility of metadynamics to sample drug access routes and diverse binding poses (Fig. 1, B–D). However, these previous studies focused on bacterial or homology-based Na_v structures and/or only examined a small number of drugs. Thus, we do not have a systematic structural comparison of the different Na_v inhibitor types within a relevant mammalian Na_v channel pore cavity.

Despite the growing number of structural and simulation studies, it is not clear if these compounds have a single stable binding site as seen in the cryo-EM structures, or multiple low-energy binding poses as characterized in previous simulations. Although there is evidence that drugs can pass the fenestrations in bacterial Na_v channels, we don’t clearly know how each inhibitor accesses the cavity binding site from the extracellular space, especially in human Na_v channels that feature four non-identical fenestrations. Furthermore, little is known about how the different molecular structures of drugs give rise to differences in pore-binding interactions or even the most likely protonation state of many of these drugs in the pore. Consequently, we require a more thorough understanding of protein–drug interactions, especially within the enigmatic ion channel pore cavity, to better manage the numerous cases of refractory channelopathies and adverse reactions to certain inhibitors.

Here, we make use of metadynamics, a biased molecular dynamics simulation method, to investigate the access routes and binding modes of nine disparate inhibitors that encompass a range of therapeutic profiles (Table 2). We examine both the neutral forms of these nine compounds, as well as the protonated forms of four of these drugs that are likely to be protonated in physiological conditions. We first show that the biased sampling method can reproduce the binding site for flecainide in the cardiac Na_v channel, as captured in the cryo-EM structure, 6uz0 (Jiang et al., 2020). Subsequently, we sample the binding of eight other inhibitors within the pore and show that compounds can differentially adopt diverse binding modes, located centrally in the pore or more peripherally in the fenestrations.

All compounds were assessed using the Major Microspecies Plugin in MarvinSketch version 23.12, 2023, ChemAxon (http://www.chemaxon.com) to determine the sites of protonation and the ratio of neutral to protonated species at pH 7. For flecainide, mexiletine, bupivacaine, and etidocaine, a restrained electrostatic potential (RESP) charge calculation was performed using Gaussian 16 Revision C.01 (RRID:SCR_014897) to obtain the charge distributions for protonated versions of these compounds. Parameterization of neutral were conducted using AnteChamber to supplement any missing definitions in the Amber GAFF2 forcefield (Wang et al., 2004, 2006).

The rat Na_v1.5 structure containing flecainide (PDB ID: 6uz0) was used to assess pore cavity accessibility and binding of different drug molecules. Voltage sensors and extraneous extracellular loop regions were removed from the structure, leaving only the transmembrane pore region (residues: 236–285, 313–417, 824–945, 1321–1504, 1641–1780). Flecainide and other co-resolved ligands were also deleted. The pore construct was embedded in a POPC membrane using PACKMOL-MEMGEN. The size of the periodic box measured 96 × 96 × 122 Å, and the system was solvated in a 0.15 M concentration of NaCl using the OPC water model (Izadi et al., 2014) and the corresponding 12-6 Li, Merz ion parameters (Li et al., 2015; Sengupta et al., 2021). Simulations were run with Assisted Model Building with Energy Refinement (AMBER2020) (RRID:SCR_014230) (Case et al., 2020) using Amber ffSB19 protein (Tian et al., 2020) and Lipid21 (Dickson et al., 2022) forcefields. The 3D molecule structures of the compounds were obtained from PubChem (RRID:SCR_004284) and manually positioned within the pore cavity space using visual molecular dynamics (VMD) (RRID:SCR_001820) (Humphrey et al., 1996), ensuring no drug–protein atomic overlaps. 20 replicate systems were each equilibrated under minimization, gradual heating, and pressurizing steps, followed by a short 5-ns simulation with the protein backbone restrained at 5 kJ/mol, and then a gradual reduction of restraint over 24 ns. 20 equilibrium simulations were run for 100 ns each using hydrogen mass repartitioning (Hopkins et al., 2015) to allow a 4-fs timestep. Simulations were maintained at 310 K and 1 bar using the Langevin thermostat (Loncharich et al., 1992) and the Monte Carlo barostat (Åqvist et al., 2004), respectively.

Well-tempered metadynamics was run using 20 walkers, with each walker initiated using the endpoints of the 20 replicates from unbiased simulation. This current protocol adapts a similar metadynamics procedure used in a previous study (Tao and Corry, 2022). PLUMED 2.7.1 (RRID:SCR_021952) (Bonomi et al., 2009; Tribello et al., 2014) was used to specify the x, y, and z coordinates of the center of mass (COM) of the drug as a 3D collective variable definition for biasing drug exploration of the protein cavity space. The drug COM position was defined relative to the central axis of the Na_v channel pore cavity using PROJECTION_ON_AXIS to determine the position along and perpendicular distance from the central pore axis and TORSION to compute the rotation around the axis. A depiction of the axis relative to the pore module and rotation angle around the axis are shown in Fig. 1, B and D. The top and bottom of the central axis were constructed by finding the COM of 100 α-carbon atoms at the extracellular and intracellular sides of the pore module, respectively. The z coordinate corresponded to the position along the axis, while x and y coordinates were calculated by multiplying the projection from the axis with the cosine and sine of the angle of rotation, respectively. The pore domain was also aligned at each timestep based on the transmembrane helix α-carbon selection.

Metadynamics parameters for biasing were applied as follows: 0.6 Å Sigma-Aldrich (or Gaussian width along x, y and z), 1.1 kJ/mol hill height, 10 bias factor (for reduction of hills), hill deposition at every 250 ps (or 1,000 steps). We found that these values were optimal for ensuring adequate sampling of the pore cavity and fenestration space within a feasible simulation timeframe.

Convergence of metadynamics simulations was assessed by comparing the FES produced at different time points. Minimal differences were seen in the total surface explored and location of minima at 2,000; 3,000; and 4,000 ns total time for neutral flecainide (Fig. S10 A). Although different drugs may vary in the time required to converge, the number of FES grid points in the pore and fenestration region sampled over the course of metadynamics simulations plateaus at around 2,000 ns (Fig. S10 B). Thus, we justified that 2,000 ns was an appropriate timeframe to produce reasonably converged surfaces across the board.

The FES for each drug was reconstructed from metadynamics hills files using the PLUMED sum_hills function. 2D surfaces for top-down and transverse views of the pore were extracted by integrating each of the third dimensions, respectively.

The top five clusters of drug binding poses were extracted using the WMC PhysBio clustering function in VMD with an RMSD cut-off of 2.5 Å across all frames of the 20 replicates of equilibrium simulation trajectories combined, striding every 0.5 ns. In the test case of flecainide, backbone residues in the pore were aligned followed by calculation of the RMSD of the heavy atoms of neutral and protonated flecainide clusters compared to the cryo-EM flecainide bound structure. MDAnalysis (RRID:SCR_025610) (Michaud-Agrawal et al., 2011) was used to compute the average structure of each cluster to find the most representative frame for each of the five clusters. COM coordinates of the five representative poses were determined using the post-processing PLUMED driver function to assess their overlap with FES minima.

The top five clusters were pooled into one trajectory and processed using ProLIF (Bouysset and Fiorucci, 2021) to identify the dominant interactions between the drug and protein residues. Default ProLIF interaction parameters were used. Van der Waals interactions were excluded. Residues were shown only if they were identified to interact with the drug in >30% of the top five clustered frames.

To assess the fenestration accessibility by different drugs, the three-dimensional FES for each drug was divided into five key regions corresponding to the central cavity and the four fenestrations, DI–II, DII–III, DIII–IV, and DI–IV. Each free energy value in these regions (0.4 Å apart) was transformed to a probability value via the equation, p = e^−ΔG/kT, where kT = 2.577 and ΔG is the free energy value. To calculate the percentage of occupancy in the four different fenestrations compared to the central cavity, we determined the sum of all the probability values in each region and divided by the total of all the regions combined.

Lipid occupancy within the fenestrations was determined by computing the volumetric density across trajectories of lipid tail atoms using cpptraj (Roe and Cheatham, 2013) volmap (at 0.4 Å resolution) and then flattening the 3D density map by averaging across the z dimension. Density values within each of the fenestration regions were averaged to compare the degree of lipid occupancy between the four fenestrations.

To estimate the heterogeneity within the cryo-EM density map of flecainide resolved in the pore of rNa_v1.5 (PDB accession no. 6uz0, EMD: accession no. 20949), OccuPy (Forsberg et al., 2023) was used to determine the local scale value of the density modeled 2.5 Å around the modeled flecainide heavy atoms. Kernel settings for OccuPy were: lowpass = 8 Å, radius = 2 pixels, size = 5 pixels, τ = 92%, tile size = 12 pixels. Scale map was visualized in UCSF ChimeraX (RRID:SCR_015872) (Meng et al., 2023).

Fig. S1 shows the RMSD of all drugs and the percentage of frames in the top five clusters. Fig. S2 and Fig. S3 show the 2D free energy surfaces for all the drugs simulated, in top-down and transverse views respectively. Fig. S4 and Fig. S5 show the representative poses of the top five clusters of the drugs in the pore, in top-down and transverse views, respectively. Fig. S6 and Fig. S7 show the residue interactions identified across the top five clusters for the different drugs. Fig. S8 shows analyses of potential correlations between drug promiscuity/% pore occupancy and physiochemical/pharmacokinetic properties. Fig. S9 shows a sequence alignment of the pore across Na_v subtypes. Fig. S10 shows the convergence of metadynamics simulations.

Metadynamics was used to bias the position of pore-blocking drugs relative to the Na_v pore module. This was done by calculating the projection along the central axis of the pore to determine the z position of the drug’s center of mass (COM) (Fig. 1 B). Polar coordinates were used to define the COM position as the drug rotates about the central axis of the channel (Fig. 1 D)—these values were then converted to x and y cartesian coordinates. Drug sampling was enhanced by biasing the x, y, and z positions of the COM.

To verify the ability of the metadynamics simulations to predict drug binding, we first tested if we could reproduce the cryo-EM–resolved binding pose of flecainide in Na_v1.5 (Jiang et al., 2020). Simulations of both neutral and protonated flecainide in the Na_v1.5 pore were able to effectively identify a binding location and poses similar to what is seen in the original cryo-EM structure (6uz0). The free energy surface (FES) derived from the simulations shows a broad free energy minimum near pore-lining S5 residues of domains II and III (Fig. 2, A and B). The location of this minima overlaps with the position of flecainide captured in the original cryo-EM structure (indicated by +). Notably, however, the simulations suggest flexibility within the binding pose as indicated by the broad minimum and multiple similarly favorable energy minima.

Clustering the simulation frames indicates the most likely binding poses for each case. This reinforces that neutral flecainide tends to bind more peripherally in the Na_v1.5 pore cavity with four of the top five clusters on the outer side of the cavity and one cluster even showing binding deep in the DII–III fenestration (Cluster 3, Fig. 2 A). Protonated flecainide shows a more central binding position in the pore (Fig. 2 B) with the representative pose from the top cluster (Fig. 2 C, red) aligning closely with the cryo-EM orientation (grey) with an RMSD of 3.2 Å (Fig. 2 D). In comparison to the neutral flecainide position, the protonation on the piperidine ring orients upwards toward the negatively charged selectivity filter as previously proposed for charged inhibitors (Buyan et al., 2018). This supports that the cryo-EM structure captures a protonated flecainide state, likely accessing the binding site via a use-dependent mechanism.

While our top binding position for protonated flecainide aligns closely with the cryo-EM structure, our simulations suggest a diversity of nearby binding poses of similar free energy. Given this, we see a greater number of interacting residues compared with that seen in the cryo-EM structure (Fig. 2 E). The binding of both versions of flecainide is stabilized by predominantly hydrophobic interactions with V933, F1420, F1461, and L1464, similar to the interactions identified in cryo-EM). Additionally, L410, F937, I1468, and F1762 also appear to be important interacting residues in our simulations. Notably, all top five clusters of neutral flecainide clusters and four out of five clusters for protonated flecainide feature the same trifluoromethylated arm in the DII–III fenestration (Fig. 2 F) aligning with the strongest density in the cryo-EM map (Fig. 2, G and H). This suggests that this arm anchors the drug while other parts of the compound can adopt multiple related poses.

We next use the metadynamics protocol to examine the binding and access of eight other pore-blocking compounds. Metadynamics effectively increased the positions around the pore explored by all drugs compared with unbiased simulations (Fig. S1 A). Given the enclosed nature of the pore and absence of deep binding cavities, the binding of inhibitor molecules is generally promiscuous and with multiple possible low-affinity binding sites for all drugs simulated (Fig. S2 and Fig. S3). Since distinct binding sites are difficult to characterize, we divided the pore into separate regions by assessing the top five clusters of each drug (Fig. S4 and Fig. S5). The majority of compounds simulated favored binding locations within the central pore cavity (Fig. S2). Some other inhibitors are shown to favor binding outside of the central pore space, with binding sites identified within the lateral fenestrations (Fig. S2). Overall, larger and protonated forms of compounds bind centrally within the pore cavity, whereas smaller hydrophobic molecules bind in the fenestrations.

To better characterize the different binding classes, we closely examined some distinct examples of pore binding (Fig. 3), mixed binding (Fig. 4), and fenestration binding (Fig. 5), as well as compared the difference between neutral and protonated versions of the drugs. Of the compounds simulated, bupivacaine displayed a strong preference for occupying the central pore and limited occupancy of the fenestrations (Fig. 3, A–D). Neutral bupivacaine has the most clearly defined central binding minimum in the FES; however, the clustering suggests it is much more flexible in this position and can orientate in a variety of ways (Fig. 3, A and C).

For mexiletine, a mix of fenestration and pore binding modes is seen. Binding pose clusters occupied an intermediary zone, whereby one or more clusters fell outside of 5 Å of the central axis with distinct minima in the peripheral pore and the fenestrations (Fig. 4, A–D). Interestingly, while both the neutral and charged versions showed a mix of pore and fenestration binding there were distinct differences between the two. The top binding pose for neutral mexiletine is seen in the DII–III fenestration (Fig. 4 C); however, this site is absent for charged mexiletine as it has a greater preference for binding in the pore or the DI–II fenestration (Fig. 4 D).

In contrast, riluzole and lamotrigine exhibit a distinct preference for fenestration binding (Fig. 5), as no COM positions of the top five clusters fall within 5 Å of the central axis of the pore. These compounds can occupy any of the fenestration except DI–IV, previously shown to be the most occluded by aromatic side chains (Tao and Corry, 2022). These compounds are the most hydrophobic and elongated and are well accommodated by the narrow fenestrations. Similar to flecainide, riluzole contains a trifluoromethyl group, which is seen to favorably occupy the DII–III fenestration.

To characterize the differences between the pore and fenestration binding modes, we determined the common residues that multiple drugs interact with. Across all the drugs studied, we observe the greatest frequency of interactions on DII and DIII (Fig. S6 A). Three hydrophobic residues lining the DII–III fenestration, V933, F937, and F1461, are shown to interact frequently with the neutral versions of flecainide, bupivacaine, and mexiletine across several of the top poses (Fig. S6 and Fig. S7). The aromatic residue F1420, located on the DIII pore loop, is also a frequently interacting partner. Given its more central location, protonated compounds interact more frequently with F1420 than their neutral counterparts. Additional interacting residues, L847, F895, and L898, are unique to neutral mexiletine due to their binding in the DI–II fenestration (Fig. 4 E and Fig. S7 B). F1762 is in the DIII–IV fenestration, and its interaction is seen for flecainide, bupivacaine, and lamotrigine, which all bind in or near the mouth of this fenestration (Fig. 2 E, Fig. 3 E, and Fig. 5 E). Overall, protonated compounds show a greater preference for interacting with polar residues on the DI pore loop, T371 and Q372 (Fig. 3 E, Fig. 4 E, Fig. S6 A, and Fig. S7).

Given its central binding position, bupivacaine is stabilized by hydrophobic interactions with residues on the pore helices of all four domains (Fig. 3 E). The majority of these interactions are hydrophobic in character, helping to explain the diversity of poses that can be adopted. The charged, protonated form of bupivacaine displays the most concentrated set of binding poses, with stronger interactions with the aromatic residues on DIV and with the charged portion angled toward the selectivity filter (Fig. S7 C). Bupivacaine shows clusters and minima lower in the cavity, with an increasing number of interactions at the cytoplasimic end of the S6 helices, particularly with residues involved in intracellular gating.

Across all drugs, their centrally located clusters are seen to interact with the hydrophobic residues at the base of the S6 helices, L410, F937, I1468, and to a lesser extent Y1769 (Fig. 2 E, Fig. 3 E, and Fig. S6 A). These residues are positioned lower in the cavity, facing toward the intracellular gate (Fig. S6 B). The presence of a drug associated with these residues could contribute to the stabilization of the closed gate.

The disparate binding of the inhibitors is predominantly stabilized via non-specific hydrophobic interactions with Na_v channel pore cavity and fenestration residues (Fig. S6 C). Although instances of other interactions were detected throughout the clusters, including pi-stacking, cation-pi, hydrogen bonding, and halogen bonding, these contribute to inhibitor binding less significantly than the hydrophobic interactions (Fig. S6 C).

To quantify the ability of different drugs to access the four fenestrations compared with the central pore, values of each drug’s 3D FES were converted to probabilities and summated over each region defined to be the pore and the four fenestrations. By comparing percentages of pore and fenestration occupancy (Fig. 6 A), we can distinguish between drugs that favor pore occupation compared with the fenestrations.

All compounds can access at least one of the four fenestrations, and most inhibitors pass through three of the fenestrations with relative ease (Fig. S2 and Fig. S3). The DII–III fenestration is the most occupied by the different drugs, followed by the DI–II fenestration in some cases (Fig. 6 A). Compounds were unable to access the DI–IV fenestration, which aligns with previous findings (Tao and Corry, 2022) that this fenestration is restricted by bulky hydrophobic aromatic residues such as F398, F402, and Y1767 in human Na_v1.5 and unfavorable for drug access. The lack of access to the DIII–-IV fenestration was unexpected, given that it was previously shown to adopt a very open conformation. We hypothesize that lipids compete with drugs in accessing the fenestration. To investigate this, we determined a density map of all the lipid tail atoms across the simulations (2D representation visualized in Fig. 6 B) and compared the density in each fenestration. We see that all four fenestrations can be occupied by lipid tails to some extent (Fig. 6, B and C); however, the darkest and most occupied region, on average, is the DIII–IV fenestration (Fig. 6, B and E), suggesting that lipids are the most stably bound in this fenestration, thus limiting the access of drugs. In the case of the DI–II and DII–III fenestrations, lipids are more mobile and can move out of these fenestrations to accommodate a drug. A representative example of this is shown in a snapshot where riluzole displaces lipid tails in the DII–III fenestration (Fig. 6 D).

Compounds of lower molecular weight, such as mexiletine, lamotrigine, lacosamide, and riluzole, sample a larger volume of pore cavity space compared with larger compounds. Given the low barriers in the FES, these drugs easily access and bind the DI–II and DII–III fenestrations. Larger compounds have greater difficulty passing through fenestrations, as seen for bupivacaine, sipatrigine, and etidocaine. Particularly, the bulky rigid compound sipatrigine has restricted movement around the pore, as suggested by the FES, and more difficulty accessing fenestrations.

For mexiletine, tolperisone, and lacosamide, with a strong indication of binding in the fenestrations, they also displayed one or more binding positions far out in the fenestrations but with a low energy barrier separating the pore region from the fenestration–membrane region. This suggests that certain compounds have a greater capacity to move in or out of the central pore region.

Although binding in the sodium channel pore is a common feature of many local anesthetic, anticonvulsant, and antiarrhythmic drugs, we still lack a detailed picture of how these drugs bind to inhibit the sodium channel. In this work, we use a molecular simulation approach to assess the dynamic binding modes of nine different pore-binding inhibitors of the sodium channel to better understand the nuances of how these pore-blocking drugs behave. We first show that our enhanced sampling metadynamics simulations can reproduce the binding site of flecainide, comparable with what is seen in cryo-EM structures. Furthermore, we confirm that the cryo-EM structure likely captures the protonated version of the drug.

Using the metadynamics protocol, we then show that different compounds have diverse binding modes within the sodium channel pore, ranging from those that occupy the central space of the pore cavity to those that prefer to bind more peripherally in the lateral fenestrations (Fig. 6, F and G). Compounds that are larger in size and protonated compounds prefer to bind in the central cavity, or have a mixed binding mode where they favor one side of the pore and bind near the mouth of a fenestration. Drugs that bind within the fenestrations are typically smaller and more hydrophobic and are stabilized by hydrophobic interactions with various fenestration-lining residues. Lipid tails also occupy the fenestrations and are seen to compete with drug binding. The DI–IV fenestration remains most restricted to drug access due to occlusion by aromatic residues. The DII–III appears to be the most commonly accessed fenestration, as various drugs can outcompete lipid tails to bind in this fenestration over the others.

The sodium channel fenestrations are seen to be occupied by lipids in molecular simulations, and lipid-like densities have been identified in cryo-EM structures. Given the hydrophobic nature of the fenestrations, this is not unexpected. However, we see that the lipids are mobile and can easily leave and re-enter some fenestrations, whereas they are more stable in other cases. Hydrophobic drug molecules are able to displace lipid occupation to access certain fenestrations. This is most apparent for the DII–III and DI–II fenestrations.

Growing evidence now suggests that the lateral fenestrations in the pore domain of ion channels are critical for their function, as well as being important regions for drug targeting (Ghovanloo et al., 2021; Huang et al., 2023a; Li et al., 2022, 2024; Hollingworth et al., 2024). Recent works hypothesize that residues in the fenestrations are involved in channel gating (Gamal El-Din and Lenaeus, 2022), and therefore if a drug can bind in the fenestrations, it could impact the ability to transition between open and closed states. Previously, we explored whether the fenestrations could provide avenues for achieving subtype-selective drug access; however, this was shown to be unlikely given their conserved nature (Tao and Corry, 2022). Nonetheless, targeting the fenestrations for drug binding could represent a novel way of inhibiting channels by carefully modulating transitions between open and closed pores without necessarily occluding sodium ion conduction.

Fenestration access is critical for the inhibitory activity of compounds like cannabidiol in Na_v1.4 (Ghovanloo et al., 2021). The authors show that cannabidiol binds near the fenestrations to stabilize the inactivated state, and mutations that restrict the fenestrations lead to decreases in cannabidiol potency. Recent experiments (Földi et al., 2021) also highlight the unique ability of riluzole to modulate pore gating but not block the conduction pathway. Fenestration-binding drugs have also been explored as a novel therapeutic approach to target loss-of-function variants of Na_v1.5 (Abdelsayed et al., 2022).

The fenestrations may serve as the binding destination of some drugs, whereas in other cases, they act as transient or metastable binding sites before the drug ultimately binds in the cavity. This bears resemblance to G protein-coupled receptors (Fronik et al., 2017), where drugs can bind at metastable secondary sites as well as the primary binding site. Given the disparity of binding locations in the Na_v channel pore, it is feasible that more than one drug molecule can concurrently bind. For example, lamotrigine may bind simultaneously in both DI–II and DII–III fenestrations and potentially lead to longer-acting inhibition. Previous work showed that Na_v channels were effectively blocked by a bivalent ligand consisting of two lidocaine molecules joined by a linker (Smith et al., 2006). However, it is yet to be explored whether two different compounds can be linked to block Na_v channels. Exploiting the multivalent interactions in the pore, by simultaneously targeting both the fenestrations and the central cavity, could be a way to enhance the specificity of Na_v channel inhibition.

Although cryo-EM structures tend to represent pore-blocking drugs adopting one distinct binding site, mutagenesis and computational studies suggest a greater diversity in binding positions. Two DIV S6 residues, F1762 and Y1769, are implicated in local anesthetic and antiarrhythmic drug binding; however, other hydrophobic residues located on DI S6 and DIII S6 have also been shown to reduce drug binding when mutated (Yarov-Yarovoy et al., 2001; Ragsdale et al., 1994, 1996; Qu et al., 1995). These data suggest that the drugs may share some common interactions but orient to different regions of the pore to prevent channel function. Unbiased flooding simulations of the human Na_v1.5 channel on the microsecond timescale further showed that neutral lidocaine has two distinct binding sites in the pore cavity (Nguyen et al., 2019). Using other enhanced simulation approaches, previous studies also show that pore inhibitors bind promiscuously, with several energy minima in the pore, rather than a single low-energy position (Buyan et al., 2018, 2021). Prior work in bacterial sodium channels also supports the idea of multiple binding sites within the pore (Boiteux et al., 2014).

Our work further supports the idea that a drug may adopt different pore-binding sites depending on whether the drug is in the neutral or protonated form. Previous work highlighted that neutral and charged compounds have distinct binding poses, with the positively charged portion of the drugs extending into the selectivity filter (Buyan et al., 2018). Our work also shows similar binding differences between neutral and protonated compounds like flecainide, where the charged region favorably orientates toward the selectivity filter and reduces the diversity of binding poses adopted.

To investigate the possible cause and effect of drug promiscuity, we calculated correlations with drug physiochemical properties and a set of consistent experimental values from Lenkey et al. (2010) (Table 2). The strongest correlation is between logP (hydrophobicity) and promiscuity (R² = 0.765, Fig. S8 B), suggesting that the more hydrophobic the drug, the lower the promiscuity of binding within the Na_v cavity. Promiscuity of the neutral compounds also appears associated with IC₅₀ (R² = 0.64, Fig. S8 D) and affinity to the inactivated state (R² = 0.652, Fig. S8 F), where decreased promiscuity is correlated with greater potency. Greater drug promiscuity is also associated with slower kinetics of drug association and dissociation, but only when we account for protonation (Fig. S8 E). This suggests that despite protonated species having greater promiscuity, they may access the pore cavity more gradually and remain confined in the cavity for longer.

Recent cryo-EM structures of pore binding in the Na_v channel cavity model a single binding site in the pore cavity for each of flecainide, riluzole, and lamotrigine (Li et al., 2024; Wu et al., 2023; Huang et al., 2023a, 2023b). This contrasts with the multiple similar sites we see in our simulations. This may be due to a range of factors affecting the interpretation of both the experimental and simulation data. First, overlapping binding sites are difficult to differentiate within the resolution of typical cryo-EM studies, for which a single model is usually fit to the density data. For flecainide, for example, the strongest density seen in the cryo-EM map is at the location of the common aspect of the poses seen in our simulations, and the more diffuse density at other regions of the compound may be caused by it occupying a diversity of positions. Lower probability sites will also produce weaker densities that are difficult to differentiate from other molecules such as water or lipids in the cavity.

To quantify the resolvability of drugs in the pore, we checked the Q-score (Pintilie et al., 2020) of drugs that have been captured in the pore of recent Na_v channel cryo-EM structures (Table 1). Comparing the ligand Q-score averages to the average Q-score of residues in the pore helices, we find that the drug is less resolved than the pore residues for the majority of structures.

Another complicating factor can be the temperatures used during cryo-EM sample preparations that may alter the relative probability of finding the drug at each position. The plunge-freezing process during cryo-EM sample preparation causes a decrease in conformational heterogeneity of biomolecules, leading to high-resolution structures but not representative of the structural ensemble at physiological temperatures (Bock and Grubmüller, 2022). For example, the energy barriers between some of the poses here are small enough that the relative population in each can be influenced by the freezing process. While the temperature dependence of conformational dynamics has been considered for other proteins and thermosensitive channels, it has not been explored in the case of drug binding in sodium channels. Therefore, it remains uncertain whether the dynamics of channel conformation and drug binding in sodium channels are affected by temperature.

On the simulation side, our results are limited by the forcefield accuracy, particularly the uncertainty of drug parameters. While we used established methods for parameterizing the drugs, it is possible that the choice of parameters or force field (e.g., if we had used CHARMM rather than AMBER force fields) could influence our results. Furthermore, since we only simulated one enantiomer for each drug, in the case of drugs with chiral centers (i.e., flecainide, mexiletine, bupivacaine, etidocaine, and tolperisone), our simulations do not account for the possibility of stereospecific interactions within the pore. Since compounds such as lacosamide (Rogawski et al., 2015), mexiletine (De Luca et al., 1995), and bupivacaine (Mitra and Chopra, 2011; Valenzuela et al., 1995) have stereoselective effects when it comes to clinical potency and safety profiles, it would be interesting to explore binding differences between enantiomers. We used metadynamics to accelerate the sampling of the pore by the drug; however, we could not guarantee whether the dynamics of the Na_v channel itself or the surrounding lipids have been sampled. Notably, our simulations started from the presumed inactivated conformation of the pore and one membrane conformation containing only POPC, therefore, it is unclear whether drug binding and accessibility deviate from other conformational states or membrane compositions. An improved structural understanding of the resting and open states of the sodium channel is warranted to better examine the molecular bases of tonic and use-dependent binding differences. Although our study only tests a specific Na_v channel subtype, rat Na_v1.5, the highly conserved nature of the pore (Fig. S9) likely means that these findings are generalizable to the other subtypes and the human Nav channels.

While we closely reproduce the cryo-EM binding pose of flecainide, the recently resolved binding sites of riluzole (PDB ID: 8thg) and lamotrigine (PDB ID: 8thh) in the Na_v1.7 pore (Huang et al., 2023b) do not match the top five clusters we identified in our metadynamics simulations (Fig. 5, C and D). Given that riluzole and lamotrigine are fenestration-occupying compounds, we found that the presence of lipids in the fenestrations competes with drug binding. That is, the nature of the lipids or detergent present will influence the most likely pose of fenestration-binding drugs. Although both riluzole and lamotrigine can favorably displace lipids in the DI–II and DII–III fenestrations, lamotrigine has a greater propensity to compete with lipid tails in the DIII–IV fenestration and is subsequently able to bind in this fenestration (Fig. 5, A–D). In fact, our data suggest that lamotrigine can bind more deeply in the DIII–IV fenestration than what is identified in the cryo-EM pose.

The disparity of chemical characteristics on the pore-blocking compounds indicates the promiscuous nature of binding in the pore, more akin to that seen in certain multidrug and peptide transporters than in classical protein-drug binding sites (Du et al., 2018). Such transporters feature large polyspecific substrate pockets, with an abundance of hydrophobic residues and occasionally a charged or polar residue, capable of accommodating diverse compounds, in a variety of orientations and positions. In the case of Na_v1.5, we found that compounds can also fluctuate in the confined cavity, finding their own favored interactions with several common residue partners, while making use of anchor points such as the charge of the selectivity filter (for protonated compounds), the aromatic cavity residues (for aryl compounds), or the hydrophobic fenestrations (for small hydrophobic compounds).

In summary, here we characterize the behavior of nine different sodium channel inhibitors within the pore of Na_v1.5, showing that there is a diversity of binding positions sampled by different drugs and that even a single drug can move between numerous similarly probable positions in the pore. While this contrasts with the traditional notions of a drug having a well-defined binding position on its receptor, this better explains the known behavior of voltage-gated sodium channel inhibitors. That is, diverse residues have been shown to be associated with drug binding, these drugs are reversible inhibitors, and above all, a great diversity of compounds are known to inhibit these proteins to exhibit a variety of effects. The subtle variations in binding help explain the disparate nature of pore inhibition and will help to guide future drug discovery efforts approaching this task from the perspective of polyspecificity rather than discrete rigid binding sites.

The data underlying free energy surfaces and interaction plots are available on Zenodo at http://doi.org/10.5281/zenodo.10408300.

Jeanne M. Nerbonne served as editor.

This research was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI) provided by the Australian Government, accessed through the national and ANU merit allocation schemes.

This work was also supported by the Australian Government Research Training Program (RTP) Scholarship.

Author contributions: E. Tao: Conceptualization, Formal analysis, Investigation, Methodology, Visualization, Writing - original draft, Writing - review & editing, B. Corry: Conceptualization, Supervision, Writing - review & editing.