Physical basis of specificity and delayed binding of a subtype selective sodium channel inhibitor

Nerve and muscle signalling is controlled by voltage-gated sodium (Nav) channels which are the targets of local anesthetics, anti-epileptics and anti-arrythmics. Current medications do not selectively target specific types of Nav found in the body, but compounds that do so have the potential to be breakthrough treatments for chronic pain, epilepsy and other neuronal disorders. We use long computer simulations totaling more than 26 μs to show how a promising lead compound can target one Nav implicated in pain perception and specific channels found in bacteria, and accurately predict the affinity of the compound to different channel types. Most importantly, we provide two explanations for the slow kinetics of this class of compound that limits their therapeutic utility. Firstly, the negative charge on the compound is essential for high affinity binding but is also responsible for energetic barriers that slow binding. Secondly, the compound has to undergo a conformational reorientation during the binding process. This knowledge aids the design of compounds affecting specific eukaryotic and bacterial channels and suggests routes for future drug development.

: Summary of the simulations run for this study.
Hunting for binding sites PFZ was placed at high concentration (69 mM) in the aqueous region of the simulation as shown in Fig. S3 and then allowed to freely interact with the protein for 1060 ns (so called 'flooding simulations'). The most commonly occupied positions, specific drug-protein interactions and free energies are characterised. To do this, the simulation system is divided into 2x2x2Å grid squares and the probability of finding the centre of mass of a PFZ molecule in each box is determined. Free energy plots were derived from this, taking advantage of the four-fold symmetry of the system, using Where P(x,y,z) is the normalised probability of finding PFZ in each grid box. two dimensional projections of the data are made by integrating the probability across one dimension between two bounds and recalculating the free energy as above.
Access to the crystallographic site The free energy of drug entry is determined using the method of umbrella sampling. [17] For the onedimension PMFs, the distance of the negatively charged warhead from the backbone of the R4 gating charge is used as the collective variable. This was defined as the distance of the centre of mass of the 9 heavy atoms of the 'warhead' end of PFZ, to the centre of mass of the R4 alpha carbon and the two alpha carbons each side of this. A harmonic potential with force constant 15 kcal/molÅ −2 is used to restrain the distance for a series of windows, spaced 0.5Å apart in the range 9.5-25Å for NavAb, and 9.5-30Å for Nav1.7/NavAb chimera which has a larger VS domain. An additional set of umbrella windows with force constant 25 kcal/molÅ −2 was used in the range 11-16Å for PFZ-Nav1.7/Navab chimera to ensure overlap between positions in neighbouring windows in the region where the forces are changing rapidly. To prevent having to sample all orthogonal coordinates in the bulk region the warhead was restrained in a cylinder of radius 10Å in the z-direction centred on an axis passing through the centre of mass of the 3 alpha carbons described above. The position data from each window is collectively analysed using the weighted histogram analysis method (WHAM) [18,19] to obtain the final potential of mean force (PMF) using the implementation of Grossfeld [20]. As each simulation contains 4 PFZ molecules (one in each protein subunit), there is the potential to gain 4 separate PMFs from each set of simulations. However, a complete PMF requires physical overlap between the positions sampled within each umbrella window. Although good overlap is found along the collective variable itself, in some cases the Cartesian coordinates of adjacent windows do not overlap. PMFs of specific compounds in which there is not overlap in both the CV and the Cartesian coordinates were excluded, leaving 3 complete PMFs for the PFZ-NavAb simulation and one for the PFZ-NavAb/Nav1.7 chimera, 4 for the PFZH-NavAb/Nav1.7 chimera and 1 for the PFZ-NavAb/Nav1.7 D1586A E1589Q mutant.
To determine if the umbrella windows have run for long enough to converge the PMF and to determine the appropriate amount of equilibration time to discard, we recalculate the PMF from different continuous subsets of the simulation data as shown in Figure S7. For the PFZ-Nav1.7/NavAb chimera, excluding more than the first 20 ns of data does not significantly influence the final result, but excluding less than this creates a noticeable effect. Thus, for the final plots we exclude the first 20 ns of each umbrella window as equilibration time. It is also evident that the PMF does not change as we extend the simulations beyond 80 ns suggesting that the final results have converged to within the uncertainties discussed below. A similar analysis showed that shorter simulations allowed for convergence for the other PMFs as indicated in Table S1.
Statistical uncertainty in the data can be computed using Monte-Carlo bootstrap analysis in which the PMF is recalculated for a number of random subsets of the data and the standard uncertainty of this determined. However, this yields negligible uncertainties which are barely visible on the PMF plots and do not represent the real uncertainty that comes from the slow timescale motions of the protein. A better estimate of the uncertainty can be gained by examining the reproducibility of the PMF itself. Thus, for the PFZ-NavAb system, the uncertainty is calculated from the range of values found in the 3 independent PMFs determined for the separate compounds in the simulations as described above. As we only have a single PMF for the PFZ-NavAb/Nav1.7 chimera we instead plot the maximum and minimum values of the PMF at each position found when recalculating the PMF from all possible windows of data ≥ 20 ns.
In order to calculate the 2D PMF as a function of distance from the binding site and dihedral angle, we conduct umbrella sampling with an additional bias potential on the phenyl-sulfonyl dihedral angle. Umbrella windows were constructed along the distance variable using a force constant of 20 kcal/molÅ −2 , spaced 0.5Åapart in the range 10-23Å, and while the dihedral windows were spaced 20 degrees apart with a force constant of 0.01kcal/mol rad −2 . Due to the large number of windows, each was run for 20ns with the first 5ns discarded as equilibration time prior to analysis. The PMF of the PFZ dihedral in bulk water was constructed by first equilibrating the compound in a 30 x 30 x 30Å water box with a single counter ion. Umbrella windows were spaced 20 degrees apart with a force constant of 0.01kcal/mol rad 2 . Runs lasted 7 ns per window with the first 2 ns discarded as equilibration prior to analysis.

Dissociation constants
The dissociation constant for PFZ in and NavAB/Nav1.7 chimera and in NavAb were determined by integrating the potential of mean force from bulk to the binding site following the method described previously [21,22] where K d is the dissociation constant, r is the radius of the constraining cylinder, N A is Avagadro's number, w(z) is the PMF, z 0 and z 1 are the values of the collective variable in the bulk region and binding site respectively, k is Boltzman's constant and T is the temperature.
The dissociation constant for PFZ in Nav1.5 was estimated by combining the value found for Nav1.7 with the free energy of the three mutations reported in Table 1 that converts the Nav1.7 binding site to be similar to Nav1.5.

Binding free energy changes from mutations
The change in binding free energy of PFZ from mutating specific residues is calculated by combining the free energies of alchemical transformations using the thermodynamic cycle below [23] (Eq. 3).
Using this, the change in binding free energy of PFZ generated by the mutation, ∆∆G, can be calculated from two separate simulations in which the side chain of a single amino acid is slowly morphed from one species to another. The free energy change is calculated with with PFZ present (∆G bound pert ) and without PFZ (∆G f ree pert ) and combined as follows to give the total change in binding free energy caused by the mutation: The free energy change of the mutation in the individual simulations is calculated using the method of free energy perturbation (FEP) [24]. To do this the mutation is broken into 40 λ windows of equal size run for 2 ns each (including 0.5 ns of equilibration). Each mutation is conducted in both directions and repeated 6 times to ensure reproduceability, yielding a total of 960 ns for each mutation. Uncertainties in the final values are determined from the standard error of the repeated simulations. Due to the strong interaction of PFZ with the gating charges, the compound remains in the binding site in all cases, avoiding some of the complications often seen in this kind of simulation.

Time scale for binding
The time constant for binding was caculated using transition state theory using a previously published approach [25,26]. In this, the rate constant is given by where ∆G b is the hight of the free energy barrier in the PMF and ω 0 is the frequency of the oscillations of the compound along the direction of the PMF defined as where z 0 is the location of the local minimum in the PMF on the outside of the energy barrier and m ef f is the effective mass of the compound calculated from its average velocity v 2 in simulations when at z 0 using the equipartition theorem m ef f v 2 = kT /2 Figure S1: Interaction energies of PFZ with nearby protein residues in the NavAb N67D mutant.

Additional Figures
Residues with strong interactions with PFZ are indicated.