Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Molecular Dynamics Study of Binding of µ-Conotoxin GIIIA to the Voltage-Gated Sodium Channel Nav1.4

Abstract

Homology models of mammalian voltage-gated sodium (NaV) channels based on the crystal structures of the bacterial counterparts are needed to interpret the functional data on sodium channels and understand how they operate. Such models would also be invaluable in structure-based design of therapeutics for diseases involving sodium channels such as chronic pain and heart diseases. Here we construct a homology model for the pore domain of the NaV1.4 channel and use the functional data for the binding of µ-conotoxin GIIIA to NaV1.4 to validate the model. The initial poses for the NaV1.4–GIIIA complex are obtained using the HADDOCK protein docking program, which are then refined in molecular dynamics simulations. The binding mode for the final complex is shown to be in broad agreement with the available mutagenesis data. The standard binding free energy, determined from the potential of mean force calculations, is also in good agreement with the experimental value. Because the pore domains of NaV1 channels are highly homologous, the model constructed for NaV1.4 will provide an excellent template for other NaV1 channels.

Citation: Mahdavi S, Kuyucak S (2014) Molecular Dynamics Study of Binding of µ-Conotoxin GIIIA to the Voltage-Gated Sodium Channel Nav1.4. PLoS ONE 9(8): e105300. https://doi.org/10.1371/journal.pone.0105300

Editor: Steven Barnes, Dalhousie University, Canada

Received: May 19, 2014; Accepted: July 22, 2014; Published: August 18, 2014

Copyright: © 2014 Mahdavi, Kuyucak. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: The authors confirm that all data underlying the findings are fully available without restriction. All data underlying the findings in our study are freely available in the supplemental file.

Funding: The authors have no support or funding to report.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Voltage-gated sodium (NaV) channels are responsible for initiation and propagation of action potential, which is essential for the activity of excitable cells such as neurons, heart and muscle cells [1]. NaV channels are involved in many physiological activities throughout the body, which also makes them potential drug targets for related disorders such as cardiac and neuropathic diseases [2]–[5]. Because of their prominence, much effort has gone into understanding the structure and function of NaV channels. For a very long time, no molecular structures were available for NaV channels, and the low resolution structures obtained from cryo-electron microscopy [6] were not very informative. In the absence of crystal structures, indirect methods such as studying the effect of mutations on ligand binding [7]–[12] have been the main source of information on NaV channels.

This situation has changed dramatically with the recent determination of the crystal structures of several bacterial NaV channels [13]–[17]. As in potassium channels, where solution of the bacterial K+ channel KcsA [18] started an intense period of examination of the structure-function relations, we expect a similar thing to happen in NaV channels. Already, there have been many computational studies of the bacterial NaV channels, investigating their ion permeation [19]–[27] and gating [28], [29] properties, as well as ligand binding [30], [31]. There are, as yet, no crystal structures for the mammalian NaV channels. Unlike in potassium channels, the tetrameric symmetry is lost when going from the bacterial to the mammalian NaV channels. Therefore, solution of the structures of the mammalian NaV channels is likely to be more difficult and may not follow soon. This leaves construction of homology models from the bacterial ones as the most viable route for making progress. Again this is not as straightforward as in potassium channels because, besides the loss of the tetrameric symmetry, the critical selectivity filter is not conserved either between the bacterial and the mammalian NaV channels. Nevertheless, the bacterial NaV channels provide a reasonable scaffold for construction of homology models for the mammalian ones, and such models can be constrained and validated using the large amount of functional data that have been accumulated over several decades [1]. Initial attempts in this regard include a docking study of tetrodotoxin and anesthetics binding to the pore domain of the NaV1.4 channel [32], and an MD study of Na+/K+ selectivity in NaV1 channels [33].

NaV channels are targets for many toxins, which bind with high affinity to various sites on the channel protein, disabling its normal function. For example, tetrodotoxin, saxitoxin, and µ-conotoxins bind to the channel vestibule and block the pore [1], [34]. This has been exploited in many experimental studies, where the known toxin structures were used to probe the pore domain of NaV channels [35], [36]. Because µ-conotoxins are the only peptide toxins that block NaV channels, there has been a lot of interest in their properties. The first µ-conotoxin to be characterised was µ-conotoxin GIIIA (µ-GIIIA), which selectively binds to the NaV1.4 channel [37], [38]. Its structure was determined from NMR by several groups [39]–[41]. A great deal of mutagenesis and functional studies have been performed on the NaV1.4–µ-GIIIA complex to determine the binding mode and identify the key residues involved in binding [42]–[55]. For this reason, NaV1.4–µ-GIIIA complex provides a unique system for testing the homology models of mammalian NaV channels that are constructed from the bacterial ones.

Here, we construct a homology model for the pore domain of the NaV1.4 channel based on the crystal structure of NaVAb. The stability of the NaV1.4 model is checked via molecular dynamics (MD) simulations. We then create a model for the NaV1.4–µ-GIIIA complex using the docking program HADDOCK [56], [57], followed by refinement in MD simulations. The proposed binding mode for the NaV1.4–µ-GIIIA complex is shown to give a satisfactory account of the available mutagenesis data. As a dynamical test of the complex model, we have determined the binding free energy of µ-GIIIA from its potential of mean force (PMF), which is again found to be in good agreement with the experimental value. The same computational methodology has previously been used to study toxin binding to potassium channels, where its ability to yield accurate protein-ligand complexes and binding free energies has been demonstrated [58]–[64]. The present study extends it to sodium channels, where much more work remains to be done.

Methods

Modeling of NaV1.4 Channel

We construct a homology model for the pore domain of NaV1.4 using the crystal structure of NaVAb (PDB ID: 3RVY) [13]. The sequence of NaV1.4 is taken from the Uniprot database (P15390). As shown in Figure 1, the sequences for the four domains of NaV1.4 are all quite different from that of NaVAb. This makes homology modeling of NaV1.4 not so straightforward as has been noted earlier [32]. Here, we have aligned the critical (DEKA) residues in the selectivity filter with the corresponding (EEEE) residues in 3RVY. A gap is placed between the E and W residues in domain II of the selectivity to align the conserved W residues. Then, we do multiple alignments between S5-pore-S6 sequence for all four domains of NaV1.4 and NaVAb using the program ClustalW [66]. The final alignments obtained for the P1-SF-P2 sequences are shown in Figure 1, which are the same as those in [15] but differ from the ones in [32] in the placement of the gap in domains I, III, and IV, which is after the W residues in the selectivity filter.

thumbnail
Download:
Figure 1. Alignment used in homology modeling of rNav1.4.

The DEKA residues in the four domains forming the selectivity filter (SF) are highlighted.

https://doi.org/10.1371/journal.pone.0105300.g001

While acceptable alignments are obtained for the pore domain, it is much harder to do the same for the S5-P1 and P2-S6 linker sequences in the turret because there are no good templates, and much less data are available on their function. Therefore, we have neglected these linker sequences in our current model of NaV1.4. While the S5-P1 linker faces the pore, it does not appear to be involved in binding of µ-GIIIA, hence our results are unlikely to be affected by its absence. A 3D model of the channel is created using Modeller [67] by threading the aligned NaV1.4 sequence for each domain on a corresponding domain of 3RVY. In order to refine the model and check its stability, we have performed MD simulations of the NaV1.4 model in a membrane environment. For this purpose, we have used the protocols established in previous MD simulations of ion channels [68], [69]. The NaV1.4 model is embedded in a lipid bilayer consisting of 153 POPE molecules in the x-y plane and solvated with a 0.1 M NaCl solution. Extra counter ions are included in the system to neutralize it where necessary. The system is then equilibrated in MD simulations in several stages. First the protein is fixed and the system is equilibrated with pressure coupling until the correct water and lipid densities are obtained. At that point, the x and y dimensions of the simulation box are fixed and pressure coupling is applied only in the z direction (the box size is about 95×96×84 Å3). In the second stage, the restraints on the protein atoms are relaxed gradually by first reducing those on the side chain atoms from kcal/mol/Å2 to 0 in 3 ns. Finally, the backbone atoms are relaxed in a similar manner. The resulting system is run for 0.1 µs to check the stability of the model. The rmsd of the backbone atoms of NaV1.4 formed a plateau after the first few ns which remained stable throughout the 0.1 µs of MD simulations with an average value of 0.6 Å, confirming its stability.

Modeling of NaV1.4–µ-GIIIA Complex

The conotoxin µ-GIIIA is a 22-residue peptide with the sequence RDCCTOOKKC KDRQCKOQRC CA with an amidated Ala at the C-terminal. It has three Arg, four Lys and two Asp residues with a total charge of +6e. The NMR structure of µ-GIIIA is shown in Figure 2A (PDB ID: 1TCJ) [41]. The side chains of the four basic residues involved in binding of µ-GIIIA to NaV1.4 are indicated explicitly in the figure. The structure of µ-GIIIA is stabilized by three disulfide bridges between C3–C15, C4–C20, and C10–C21. However, as seen from the superposition of the NMR snapshots in Figure 2B, the peptide has quite a bit of flexibility around the hinges at the C3 and C10 residues. For this reason, the residues 11–22 forming the binding interface will be used in the rmsd calculations when we check for distortion of the toxin during umbrella sampling MD simulations.

thumbnail
Download:
Figure 2. NMR structures of µ-GIIIA.

(A) Structure of µ-GIIIA with the pore blocking residues K11, R13 and K16 pointing downward. Three disulfide bridges and the C3–K10 hydrogen bond stabilizing the structure are indicated explicitly. (B) Superposition of the ten NMR structures demonstrating the flexibility of the N-terminal residues 1–9 around the C3 and C10 hinges.

https://doi.org/10.1371/journal.pone.0105300.g002

The initial poses for the NaV1.4–µ-GIIIA complex are generated by docking the µ-GIIIA structure to the NaV1.4 model using the program HADDOCK [56], [57]. In previous studies of toxin binding to potassium channels, HADDOCK has given very good results, reducing the time needed for refinement with MD [59]–[63], [70]. In order to get an adequate sampling of the side chain orientations, we use all ten NMR conformers of µ-GIIIA in ensemble docking. Because there are no well-known binding motifs for NaV1 channel blockers—like the pore inserting Lys in potassium channel blockers—we have considered several possibilities for restraints in HADDOCK. To facilitate comparisons with the mutation data and simplify interpretation of the results, we use a single restraint in each docking study. The EEDD and DEKA ring of residues are the potential sites on the channel for using restraints. However, the mutation data indicates that the EEDD residues play a much more important role in binding of µ-GIIIA than the DEKA residues [47]. Therefore, only the EEDD ring is used as a restraint site in the following docking studies. The potential restraint sites on the toxin include the residues R1, K11, R13, K16, and R19, which are identified in mutagenesis experiments (see Table 1). Separate docking studies are performed for each of these residues and the EEDD ring as a restraint. One thousand conformations are generated from each docking, and the top hundred selected on the basis of scoring are analyzed via clustering to find the dominant binding mode. In choosing among the five sets of complexes generated, we have also used the fact that µ-GIIIA blocks the channel. This criterion is satisfied only for the complexes in the set with the R13–EEDD restraint, which also have the highest scores. The complex structure in the top-ten of this set that has the most toxin–channel contacts is selected for refinement in MD simulations (see Table 2).

thumbnail
Download:
Table 1. Effect of the mutations in µ-GIIIA on the IC50 values for binding to Nav1.4.

https://doi.org/10.1371/journal.pone.0105300.t001

thumbnail
Download:
Table 2. List of interacting residues in the Nav1.4–µ-GIIIA complex.

https://doi.org/10.1371/journal.pone.0105300.t002

The complex structure obtained from docking is aligned with the channel model in the membrane, and the coordinates of the toxin are transferred to the channel model. The protocol used in equilibrating the channel protein is also used for the complex, with the channel protein and toxin being relaxed simultaneously. The equilibrated system is run for another 50 ns to check its stability and generate trajectory data for analysis. During this MD simulation, a small restraint with kcal/mol/Å2 is applied to channel backbone atoms to preserve the integrity of the channel but no restraints are imposed on the toxin. The system is found to be well-equilibrated from the start, and therefore, the trajectory data are used for the analysis of the complex structure.

MD Simulations and PMF Calculations

All MD simulations are performed using version 2.7 of NAMD [71] with the CHARMM36 force field [72]. An NpT ensemble is used with periodic boundary conditions. Pressure is kept at 1 atm and temperature at 300 K using Langevin coupling with damping coefficients of 5 ps−1. Lennard-Jones interactions are switched off smoothly within distance of 10–13.5 Å. Electrostatic interactions are computed without truncation using the particle-mesh Ewald algorithm. A time step of 2 fs is employed in all MD simulations. The trajectory data is saved at 1 ps intervals but the reaction coordinate is written at every time step in umbrella sampling simulations.

The PMF for dissociation of µ-GIIIA from NaV1.4 is constructed using umbrella sampling MD simulations. As the method has been described in detail previously [58], [63], we give only a brief description here. The reaction coordinate is chosen as the distance between the center of mass (COM) of the channel protein and the COM of the toxin along the channel axis. Initially 28 umbrella windows are created along the channel axis at 0.5 Å intervals using steered MD with a force constant kcal/mol/Å2 and pulling speed Å/ns. Six more windows are added subsequently to ensure that the toxin has reached the bulk region signalled by a flat PMF. The same force constant of kcal/mol/Å2 is used in umbrella sampling simulations, which is found to be optimal for toxins of this size [58]. The overlaps of distributions between the neighboring windows should be to avoid numerical instabilities in construction of the PMF. We have included two more windows between the windows 4–5 and 22–23, where this condition is not satisfied. The reaction coordinates collected from the simulations are unbiased and combined using the weighted histogram analysis method [73]. Umbrella sampling simulations are continued until the convergence of the PMF is assured from block data analysis of the PMF data.

The binding constant is determined by integrating the PMF, , along the axis [58], [63](1)where the and are the initial and final points in the PMF, and is the average cross sectional area of the binding pocket as explored by the COM of the toxin, which is determined from its transverse fluctuations. The value of is obtained from restraint-free MD simulations of the NaV1.4–µ-GIIIA complex as 0.79 Å. Finally, the standard binding free energy of µ-GIIIA is obtained from the binding constant using(2)where is the standard concentration of 1 M.

Results and Discussion

Binding Mode of the NaV1.4–µ-GIIIA Complex

Snapshots of the NaV1.4–µ-GIIIA complex obtained from docking and MD simulations are shown in Figure 3 (A PDB file giving the coordinates of the complex structure is provided in File S2). The five basic residues on the toxin that make contacts with the acidic residues on the channel are indicated explicitly. To provide a more quantitative picture of the binding mode, we have calculated the average N–O distances between the interacting residues from the 50 ns MD trajectory data (Table 2). The distances obtained from HADDOCK are also shown for comparison. Of the nine contact pairs listed, HADDOCK has got six of them within 0.3 Å of the MD result and one of them within 1.5 Å. Only in two cases, the discrepancy is larger than 5 Å. Considering that the binding mode is rather complex involving multiple partners for some of the residues, this is an excellent result from HADDOCK. It also explains why the complex structure obtained from HADDOCK has equilibrated so quickly in MD simulations.

thumbnail
Download:
Figure 3. Binding mode of the Nav1.4–µ-GIIIA complex.

For clarity, domains I and III (top) and II and IV (bottom) are shown separately. All the important interactions between the channel and toxin residues are indicated explicitly.

https://doi.org/10.1371/journal.pone.0105300.g003

In order to validate the NaV1.4–µ-GIIIA complex, we compare the binding mode characterized in Figure 3 and Table 2 with the mutagenesis data summarized in Table 1. The toxin residue R13 makes the most contacts with the channel residues and hence is expected to play a major role in binding, which is in good agreement with experiments (Table 1). The R13A mutation causes more than 200-fold increase in the IC50 ratio. From Eq. 2, this corresponds to more than 3 kcal/mol change in the binding free energy. In previous studies of binding of ShK toxin to KV1 channels, [59]–[61], we have shown that neutralizing a charged toxin residue in contact with channel residues costs about 2 kcal/mol in binding free energy, provided the binding mode is preserved after the mutation. However, if the binding mode of the mutated toxin is substantially different from that of the wild type, there could be larger changes in the binding free energy. We have, therefore, performed docking calculations for µ-GIIIA[R13A] and found that the binding mode has completely changed with R1 inserting in the NaV1.4 pore. This highlights the importance of checking the binding mode when interpreting unusual results in alanine scanning experiments. Other mutations of R13 also reduce the affinity, including the conservative mutation R13K, presumably because a Lys residue cannot sustain the multiple contacts R13 makes. R13 is followed in importance by the residues R19, K16, and K11, each of which makes tight contacts with the acidic residues on the channel. For these residues, the observed loss of affinity, when they are mutated to Ala, are mostly in the range expected from switching off the charge (Table 1). The quality of the charge contact between K8 and D1248 is not as good as the others because K8 is on the flexible N-terminal part of the toxin, which fluctuates around the C10 hinge. As a result, the K8–D1248 average distance and its sigma are substantially larger compared to the other charge contacts. Again this is consistent with the experimental observation that K8A mutation causes less affinity loss relative to the other contacts listed in Table 2 (Table 1). Further evidence for the relative strength of the individual interactions will be presented when we discuss the persistence of contacts during dissociation of µ-GIIIA from NaV1.4.

The only residue whose mutation to Ala reduces the affinity of µ-GIIIA but does not appear in the proposed binding mode is R1. As seen from Figure 3, R1 is on the opposite side of the binding interface and does not make any contacts with the channel residues (the N–O distance for the closest acidic residue on the channel is more than 9 Å). We have seen in studies of potassium channel toxins that mutation of an Arg residue could still result in a different binding mode even though it is not directly involved in binding [59], [61]. To see if a similar thing happens in µ-GIIIA, we have docked µ-GIIIA[R1A] to NaV1.4. The resulting binding mode is found to be substantially different from that of µ-GIIIA, which indicates that the reduction in affinity is likely to be caused by a change in the binding mode rather than loss of a charge interaction with the channel residues.

The binding mode emerging from this study provides a complete block of the fairly large vestibule of NaV1.4. This is illustrated in Figure 4, which shows a bottom-view of the complex. The three basic residues, R13, K16, and K11 are seen to weave around the EEDD ring, making multiple contacts with residues in all four domains. We note that there is some redundancy here because two residues such as R13 and K16 can still cover all four domains consistent with the observation that µ-GIIIA[K11A] can also block the channel [53]. The fact that two basic residues are sufficient to block the pore is also supported by several other µ-conotoxin blockers of NaV1 channels, which have only two basic residues available at the binding interface (e.g., BuIIIB and SIIIA).

thumbnail
Download:
Figure 4. Bottom view of the complex structure demonstrating blocking of the pore.

The toxin residues R13, K11, and K16 make contacts with the channel residues EEDD in all four domains.

https://doi.org/10.1371/journal.pone.0105300.g004

Compared to potassium channels there are relatively fewer blockers of sodium channels, which is due to the larger pore size in the latter. A pore inserting Lys is sufficient to block a potassium channel whereas at least two basic residues are required to achieve the same in a sodium channel. The pore inserting Lys motif has been instrumental in functional studies of potassium channels using toxins peptides as probes. This has simplified interpretation of experimental results as well as construction of computational models of channel–toxin complexes. The situation in NaV1 channels is much more complicated due to many possible configurations for coupling of 2–3 basic toxin residues with the EEDD residues in the pore. Also there are many high-affinity toxins that do not block NaV1 channels. These features have certainly made interpretation of mutation experiments a more difficult task and sometimes resulted in conflicting proposals for the binding modes. Construction of accurate complex models using homology models of NaV1 channels is expected to ameliorate this situation. Moreover such complex models will be very useful in designing analogues with enhanced affinity and selectivity properties, which may be required for development of toxin blockers of NaV1 channels as therapeutic agents.

Umbrella Sampling Simulations and Binding Free Energy

We have previously shown in over a dozen case studies involving potassium channel toxins that the binding free energy can be determined near chemical accuracy from PMF calculations [58]–[63]. Thus calculation of the standard binding free energy of µ-GIIIA will provide a complementary test for the accuracy of the proposed NaV1.4–µ-GIIIA model. The PMF for the dissociation of µ-GIIIA from NaV1.4 is constructed from umbrella sampling MD simulation as described in Methods. Distortion of a ligand during umbrella sampling simulations is a concern in PMF calculations, especially for flexible peptides, which may lead to erroneous results if the distortion becomes permanent after the ligand [65]. We have checked against this possibility by calculating the average rmsd of µ-GIIIA in each umbrella sampling window. As shown in Figure 5, the toxin undergoes some distortion while it is pulled out of the binding pocket but the elastic energy associated with this distortion is recovered once the toxin is in the bulk region as indicated by the return of the rmsd to the bulk value. Lack of distortion is also verified by aligning the µ-GIIIA structures from the last umbrella window with the NMR structure.

thumbnail
Download:
Figure 5. Average rmsd of the µ-GIIIA backbone atoms for the residues 11–22 calculated for each umbrella window.

The bulk value obtained from MD simulations of µ-GIIIA in a box of water is indicated by the dashed line.

https://doi.org/10.1371/journal.pone.0105300.g005

Two main questions in PMF calculations are how far the PMF should be extended and how long each window should be run. The first is addressed by appearance of a flat region in the PMF which indicates that the ligand has reached the bulk region. The second question is rarely addressed in PMF calculations but it is equally important because without sufficient data, one is likely to obtain a wrong answer. We address this issue by performing block data analysis of the PMF data. That is, we construct PMFs from 2 ns blocks of data and slide the blocks in 1 ns steps over the range of the available data. As shown in Figure 6, the PMFs initially drop monotonically and then fluctuate around a base line. During the first phase (1–4 ns), the PMFs drop because of the improved screening of the channel–toxin interactions as the system equilibrates. In the second phase (4–10 ns), fluctuations of the PMFs around a base line are of statistical nature and indicate that the system has been equilibrated. A common practice in PMF calculations is to exclude an arbitrary amount of data for equilibration (e.g., 1 ns), and consider the rest of the data in production of the PMF. As seen from the convergence study in Figure 6, this could result in mixing of the equilibration and production data, which would lead to an overestimation of the binding affinity.

thumbnail
Download:
Figure 6. Convergence study of the Nav1.4–µ-GIIIA PMF from block data analysis.

To minimize fluctuations, a relatively large sampling size of 2 ns is used, which is slided in 1 ns steps over the 10 ns of data. The block-data PMFs drop monotonically in the first 4 ns as the system equilibrates and then fluctuate around a base line from 4–10 ns, signalling equilibration. The final PMF obtained from 4–10 ns is indicated by a thick black line.

https://doi.org/10.1371/journal.pone.0105300.g006

The final PMF for the NaV1.4–µ-GIIIA complex is indicated by a thick black line in Figure 6. We will comment on distinct features of the PMF when we discuss the evolution of the distances between the contact pairs below. Using Eq. 1, we numerically integrate the PMF to find the binding constant and then determine the standard binding free energy using Eq. 2. The calculated value, kcal/mol, is in good agreement with the experimental value of kcal/mol, which is determined from the most recent IC50 measurement (19 nM) where state-of-the-art equipment was employed [74]. The level of agreement obtained for the standard binding free energy is similar to those obtained for potassium channel toxins [58]–[64], and provides further support for the accuracy of the proposed model of the NaV1.4–µ-GIIIA complex.

The binding mode of the NaV1.4–µ-GIIIA complex is dominated by charge interactions where the pairs are mostly at contact distances (Table 2). Analysis of how the contact distances change during dissociation provides complementary information on the relative strength of the various interactions as well as helping to explain specific features of the PMF. For this purpose, we have calculated the average pair distances in each umbrella window and plotted them as a function of the window position (Figure 7). The toxin residues are seen to break up in the following order and at the approximate positions: K11 at 30 Å, K16 at 34 Å, R19 at 35 Å, and finally R13 at 36 Å. Because the weakest interactions will break up first and the strongest ones last, the persistence length of a pair during dissociation provides a good measure for its relative strength. This intuitive picture for the strength of the interactions is in good agreement with the mutation data in Table 2. For example, according to the set of data from [53], the affinity loss for mutations of K11, K16, R19, and R13 to Ala are, respectively, 10, 81, 199, and 642.

thumbnail
Download:
Figure 7. Evolution of the distances between the interacting pairs as µ-GIIIA dissociates from Nav1.4.

Average N–O distances obtained from each umbrella window are plotted as a function of the window position. All the pairs in Table 2 are considered except K8–D1248 which dissociates immediately.

https://doi.org/10.1371/journal.pone.0105300.g007

Apart from the glitch around 31–32 Å, the final PMF in Fig 6 follows a steadily rising trajectory until all the contact pairs are broken off at 36 Å. The glitch appears to be associated with the fluctuations of the K11 and K16 pair distances (Figure 7). After Å, the charged pairs interact via screened Coulomb interactions which corresponds to the shoulder region in the PMF. For Å, the distances between the charged pairs are larger than 15 Å. At those distances, the Coulomb interactions are mostly screened off, which correlates well with the PMF leveling off after Å and becoming flat for Å (Figure 6).

Conclusions

In this work, we have constructed a homology model for the pore domain of the NaV1.4 channel by aligning the DEKA residues in the selectivity filter with the corresponding EEEE residues in NaVAb, whose crystal structure was determined recently. In order to validate the NaV1.4 model, we have used the extensive functional data obtained from binding studies of µ-conotoxin GIIIA. An initial model for the NaV1.4–µ-GIIIA complex is created using HADDOCK, which is refined in MD simulations. The binding mode obtained is in broad agreement with the available mutagenesis data and shows that the toxin blocks the pore through multiple interactions of the R13, K16 and K11 residues with the outer ring of EEDD residues in the channel vestibule. The standard binding free energy of µ-GIIIA is determined from the PMF calculations and found to agree with the experimental value within chemical accuracy. Thus the proposed model of the NaV1.4–µ-GIIIA complex has been well validated. Because there is a high degree of homology among the NaV1 channels, the present NaV1.4 model can be used as a template in constructing homology models for the pore domain of other NaV1 channels.

Our focus in this first study was the validation of the pore domain using the data on binding of µ-GIIIA. For further studies of toxin binding to NaV1 channels one needs to include the selectivity filter and the S5-P1 linker in the model. For example, to investigate the ionic strength dependence of toxin binding [54], a validated model of the selectivity filter is required. This can be achieved by studying the permeation and selectivity properties of Na+ ions, which we hope to tackle in a forthcoming paper. Our attempts to model the S5-P1 linker in the turret of the NaV1.4 channel have not been successful due to lack of good templates. Because binding of µ-GIIIA does not appear to involve the S5-P1 linker residues, a satisfactory binding mode could still be obtained without modeling this region. However, the S5-P1 linker residues are involved in binding of some other µ-conotoxins, and to understand the differences in their affinity and selectivity properties [74], it will be important to construct models of NaV1 channels including the full turret region. The available mutagenesis data could provide valuable guidance in this endeavor. Finally, there are ongoing efforts to harness the therapeutic potential of µ-conotoxins by designing analogues that selectively bind to a target NaV1 channel [3]–[5]. Construction of accurate models of NaV1–µ-conotoxin complexes will be very useful in such efforts.

After the completion of this work, a paper on binding of µ-conotoxin PIIIA to the NaV1.4 channel has appeared [75]. In this paper, two very different binding modes were predicted with almost equal binding free energies. This is very different from the unique binding mode found for µ-GIIIA in our study and needs to be investigated further. A snapshot comparing the NaV1.4 models used in the two studies is given in Figure S1 in File S1.

Supporting Information

File S1.

Figure S1, Snapshot comparing the pore domain of our Nav1.4 model with that of Chen et al. [75].

https://doi.org/10.1371/journal.pone.0105300.s001

(PDF)

File S2.

PDB file giving the coordinates of the Nav1.4-GIIIA complex model used in drawing Figure 3.

https://doi.org/10.1371/journal.pone.0105300.s002

(PDB)

Acknowledgments

Calculations were performed using the HPC facilities at the National Computational Infrastructure (Canberra).

Author Contributions

Conceived and designed the experiments: SM SK. Performed the experiments: SM. Analyzed the data: SM SK. Contributed to the writing of the manuscript: SM SK.

References

  1. 1. Hille B (2001) Ionic Channels of Excitable Membranes. 3rd ed, Sunderland, MA: Sinauer Associates.
    • 2. Ashcroft FM (2000) Ion Channels and Disease: Channelopathies. San Diego: Academic Press.
      • 3. French RJ, Terlau H (2004) Sodium channel toxins–receptor targeting and therapeutic potential. Curr Med Chem 11: 3053–3064.
      • 4. Twede VD, Miljanich G, Olivera BM, Bulaj G (2009) Neuroprotective and cardioprotective conopeptides: an emerging class of drug leads. Curr Opin Drug Discov Devel 12: 231–239.
      • 5. Norton RS (2010) µ-Conotoxins as leads in the development of new analgesics. Molecules 15: 2825–2844.
      • 6. Sato C, Ueno Y, Asai K, Takahashi K, Sato M, et al. (2001) The voltage-sensitive sodium channel is a bell-shaped molecule with several cavities. Nature 409: 1047–1051.
      • 7. Stuhmer W, Conti F, Suzuki H, Wang X, Noda M, et al. (1989) Structural parts involved in activation and inactivation of the sodium channel. Nature 339: 597–603.
      • 8. Terlau H, Heinemann SH, Sthmer W, Pusch M, Conti F, et al. (1991) Mapping the site of block by tetrodotoxin and saxitoxin of sodium channel II. FEBS Lett 293: 93–96.
      • 9. Backx PH, Yue DT, Lawrence JH, Marban E, Tomaselli GF (1992) Molecular localization of an ion-binding site within the pore of mammalian sodium-channels. Science 257: 248–251.
      • 10. Perez-Garcia MT, Chiamvimonvat N, Marban E, Tomaselli GF (1996) Structure of the sodium channel pore revealed by serial cysteine mutagenesis. Proc Natl Acad Sci USA 93: 300–304.
      • 11. Catterall WA (2000) From ionic currents to molecular mechanisms: the structure and function of voltage-gated sodium channels. Neuron 26: 13–25.
      • 12. Fozzard HA, Lipkind GM (2010) The tetrodotoxin binding site is within the outer vestibule of the sodium channel. Mar Drugs 8: 219–234.
      • 13. Payandeh J, Scheuer T, Zheng N, Catterall WA (2011) The crystal structure of a voltage-gated sodium channel. Nature 475: 353–358.
      • 14. Payandeh J, Gamal El-Din TM, Scheuer T, Zheng N, Catterall WA (2012) Crystal structure of a voltage-gated sodium channel in two potentially inactivated states. Nature 486: 135–139.
      • 15. Zhang X, Ren W, DeCaen P, Yan C, Tao X, et al. (2012) Crystal structure of an orthologue of the nachbac voltage-gated sodium channel. Nature 486: 130–134.
      • 16. McCusker EC, Bagneris C, Naylor CE, Cole AR, D'Avanzo N, et al. (2012) Structure of a bacterial voltage-gated sodium channel pore reveals mechanisms of opening and closing. Nat Commun 3: 1102.
      • 17. Shaya D, Findeisen F, Abderemane-Ali F, Arrigoni C, Wong S, et al. (2014) Structure of a prokaryotic sodium channel pore reveals essential gating elements and an outer ion binding site common to eukaryotic channels. J Mol Biol 426: 467–483.
      • 18. Doyle DA, Cabral JM, Pfuetzner RA, Kuo A, Gulbis JM, et al. (1998) The structure of the potassium channel: molecular basis of K+ conduction and selectivity. Science 280: 69–77.
      • 19. Carnevale V, Treptow W, Klein ML (2011) Sodium ion binding sites and hydration in the lumen of a bacterial ion channel from molecular dynamics simulations. J Phys Chem Lett 2: 2504–2508.
      • 20. Corry B, Thomas M (2012) Mechanism of ion permeation and selectivity in a voltage gated sodium channel. J Am Chem Soc 134: 1840–1846.
      • 21. Furini S, Domene C (2012) On conduction in a bacterial sodium channel. PLoS Comput Biol 8: e1002476.
      • 22. Qiu H, Shen R (2012) Guo W (2012) Ion solvation and structural stability in a sodium channel investigated by molecular dynamics calculations. Biochim Biophys Acta 1818: 2529–2535.
      • 23. Ke S, Zangerl EM, Stary-Weinzinger A (2013) Distinct interactions of Na+ and Ca2+ ions with the selectivity filter of the bacterial sodium channel NaVAb. Biochem Biophys Res Commun 430: 1272–1276.
      • 24. Zhang X, Xia M, Li Y, Liu H, Jiang X, et al. (2013) Analysis of the selectivity filter of the voltage-gated sodium channel NaVRh. Cell Res 23: 409–422.
      • 25. Stock L, Delemotte L, Carnevale V, Treptow W, Klein ML (2013) Conduction in a biological sodium selective channel. J Phys Chem B 117: 3782–3789.
      • 26. Ulmschneider MB, Bagneris C, McCusker EC, DeCaen PG, Delling M, et al. (2013) Molecular dynamics of ion transport through the open conformation of a bacterial voltage-gated sodium channel. Proc Natl Acad Sci USA 110: 6364–6369.
      • 27. Chakrabarti N, Ing C, Payandeh J, Zheng N, Catterall WA, et al. (2013) Catalysis of Na+ permeation in the bacterial sodium channel NaVAb. Proc Natl Acad Sci USA 110: 11331–11336.
      • 28. Yarov-Yarovoy V, Caen PG, Westenbroek RE, Pan CY, Scheuer T, et al. (2012) Structural basis for gating charge movement in the voltage sensor of a sodium channel. Proc Natl Acad Sci USA 109: E93–E102.
      • 29. Amarala C, Carnevale V, Klein ML, Treptow W (2012) Exploring conformational states of the bacterial voltage-gated sodium channel Navab via molecular dynamics simulations. Proc Natl Acad Sci USA 109: 21336–21341.
      • 30. Chen R, Chung SH (2012) Binding modes of µ-Conotoxin to the bacterial sodium channel (NaVAb). Biophys J 102: 483–488.
      • 31. Raju SG, Barber AF, LeBard DN, Klein ML, Carnevale V (2013) Exploring volatile general anesthetic binding to a closed membrane-bound bacterial voltage-gated sodium channel via computation. PLoS Comput Biol 9: e1003090.
      • 32. Tikhonov DB, Zhorov BS (2012) Architecture and pore block of eukaryotic voltage-gated sodium channels in view of NavAb bacterial sodium channel structure. Mol Pharmacol 82: 97–104.
      • 33. Xia M, Liu H, Li Y, Yan N, Gong H (2013) The mechanism of Na+/K+ selectivity in mammalian voltage-gated sodium channels based on molecular dynamics simulation. Biophys J 104: 2401–2409.
      • 34. Terlau H, Olivera B (2004) Conus venoms: a rich source of novel ion channel-targeted peptides. Physiol Rev 84: 41–68.
      • 35. Li RA, Tomaselli GF (2004) Using the deadly µ-conotoxins as probes of voltage-gated sodium channels. Toxicon 44: 117–122.
      • 36. Dutertre S, Lewis RJ (2010) Use of venom peptides to probe ion channel structure and function. J Biol Chem 285: 13315–13320.
      • 37. Cruz LJ, Gray WR, Olivera BM, Zeikus RD, Kerr L, et al. (1985) Conus Geographus toxins that discriminate between neuronal and muscle sodium channels. J Biol Chem 260: 9280–9288.
      • 38. Ohizumi Y, Minoshima S, Takahashi M, Kajiwara A, Nakamura H, et al. (1986) Geographutoxin II, a novel peptide inhibitor of na channels of skeletal muscles and autonomic nerves. J Pharmacol Exp Ther 239: 243–248.
      • 39. Ott KH, Becker S, Gordon RD, Ruterjans H (1991) Solution structure of mu-conotoxin GIIIA analyzed by 2D-NMR and distance geometry calculations. FEBS Lett 278: 160–166.
      • 40. Lancelin JM, Kohda D, Tate S, Yanagawa Y, Abe T, et al. (1991) Tertiary structure of conotoxin GIIIA in aqueous-solution. Biochemistry 30: 6908–6916.
      • 41. Wakamatsu K, Kohda D, Hatanaka H, Lancelin JM, Ishida Y, et al. (1992) Structure-activity relationships of µ-Conotoxin GIIIA: structure determination of active and inactive sodium channel blocker peptides by NMR and simulated annealing calculations. Biochemistry 31: 12577–12584.
      • 42. Sato K, Ishida Y, Wakamatsu K, Kato R, Honda H, et al. (1991) Active site of µ-Conotoxin GIIIA a peptide blocker of muscle sodium channels. J Biol Chem 266: 16989–16991.
      • 43. Becker S, Prusak-Sochaczewski E, Zamponi G, Beck-Sickinger AG, Gordon RD, et al. (1992) Action of derivatives of µ-Conotoxin GIIIA on sodium channels. Biochemistry 31: 8229–8238.
      • 44. Dudley SC, Todt H, Lipkind G, Fozzard HA (1995) A µ-Conotoxin-insensitive Na+ channel mutant: possible localization of a binding site at the outer vestibule. Biophys J 69: 1657–1665.
      • 45. Li RA, Tsushima RG, Kallen RG, Backx PH (1997) Pore Residues critical for µ-CTX binding to rat skeletal muscle Na+ channels revealed by cysteine mutagenesis. Biophys J 73: 1874–1884.
      • 46. Chang NS, French RJ, Lipkind GM, Fozzard HA, Dudley S Jr (1998) Predominant interactions between µ-Conotoxin Arg-13 and the skeletal muscle Na+ channel localized by mutant cycle analysis. Biochemistry 37: 4407–4419.
      • 47. Chahine M, Sirois J, Marcotte P, Chen LQ, Kallen RG (1998) Extrapore residues of the S5–S6 loop of domain 2 of the voltage-gated skeletal muscle sodium channel (rSkM1) contribute to the µ-Conotoxin GIIIA binding site. Biophys J 75: 236–246.
      • 48. Dudley SC Jr, Chang N, Hall J, Lipkind G, Fozzard HA, et al. (2000) µ-Conotoxin GIIIA interactions with the voltage-gated Na+ channel predict a clockwise arrangement of the domains. J Gen Physiol 116: 679–689.
      • 49. Li RA, Ennis IL, French RJ, Dudley SC Jr, Tomaselli GF, et al. (2001) Clockwise domain arrangement of the sodium channel revealed by µ-Conotoxin (GIIIA) docking orientation. J Biol Chem 276: 11072–11077.
      • 50. Li RA, Sato K, Kodama K, Kohno T, Xue T, et al. (2002) Charge conversion enables quantification of the proximity between a normally-neutral µ-Conotoxin (GIIIA) site and the Na+ channel pore. FEBS Lett 511: 159–164.
      • 51. Cummins TR, Aglieco F, Dibhajj SD (2002) Critical Molecular determinants of voltage-gated sodium channel sensitivity to µ-Conotoxins GIIIA/B. Mol Pharmacol 61: 1192–1201.
      • 52. Hui K, Lipkind G, Fozzard HA, French RJ (2002) Electrostatic and steric contributions to block of the skeletal muscle sodium channel by µ-Conotoxin. J Gen Physiol 119: 45–54.
      • 53. Xue T, Ennis IL, Sato K, French RJ, Li RA (2003) Novel interactions identified between µ-Conotoxin and the Na+ channel domain I P-loop: implications for toxin-pore binding geometry. Biophys J 85: 2299–2310.
      • 54. Li RA, Hui K, French RJ, Sato K, Henrikson CA, et al. (2003) Dependence of µ-Conotoxin block of sodium channels on ionic strength but not on the permeating [Na+]. J Biol Chem 278: 30912–30919.
      • 55. Choudhary G, Aliste MP, Tieleman DP, French RJ, Dudley SC Jr (2007) Docking of µ-Conotoxin GIIIA in the sodium channel outer vestibule. Channels 1: 344–352.
      • 56. Dominguez C, Boelens R, Bonvin AMJJ (2003) HADDOCK: a protein-protein docking approach based on biochemical or biophysical information. J Am Chem Soc 125: 1731–1737.
      • 57. de Vries SJ, van Dijk AD, Krzeminski M, van Dijk M, Thureau A, et al. (2007) HADDOCK versus HADDOCK: new features and performance of haddock 2.0 on the capri targets. Proteins 69: 726–733.
      • 58. Chen PC, Kuyucak S (2011) Accurate determination of the binding free energy for KcsA-Charybdotoxin complex from the potential of mean force calculations with restraints. Biophys J 100: 2466–2474.
      • 59. Rashid MH, Kuyucak S (2012) Affinity and selectivity of ShK toxin for the Kv1 potassium channels from free energy simulations. J Phys Chem B 116: 4812–4822.
      • 60. Pennington MW, Rashid MH, Tajhya RB, Beeton C, Kuyucak S, et al. (2012) A C-Terminally amidated analogue of ShK is a potent and selective blocker of the voltage-gated potassium channel Kv1.3. FEBS Lett 586: 3996–4001.
      • 61. Rashid MH, Heinzelmann G, Huq R, Tajhya RB, Chang SC, et al. (2013) A potent and selective peptide blocker of the Kv1.3 channel: prediction from free-energy simulations and experimental confirmation. PLoS One 8: e78712.
      • 62. Mahdavi S, Kuyucak S (2013) Why drosophila shaker K+ channel is not a good model for ligand binding to voltage-gated Kv1 channels. Biochemistry 52: 1631–1640.
      • 63. Rashid MH, Kuyucak S (2014) Free energy simulations of binding of HsTx1 toxin to Kv1 potassium channels: the basis of Kv1.3/Kv1.1 selectivity. J Phys Chem B 118: 707–716.
      • 64. Rashid MH, Huq R, Tanner MR, Chhabra S, Khoo KK, et al. (2014) A potent and Kv1.3-selective analogue of the scorpion toxin HsTX1 as a potential therapeutic for autoimmune diseases. Sci Rep 4: 4509.
      • 65. Chen PC, Kuyucak S (2009) Mechanism and energetics of charybdotoxin unbinding from a potassium channel from molecular dynamics simulations. Biophys J 96: 2577–2588.
      • 66. Thompson JD, Higgins DG, Gibson TJ (1994) CLUSTAL-W - improving the sensitivity of progressive multiple sequence alignment through sequence weighting position-specific gap penalties and weight matrix choice. Nucl Acids Res 22: 4673–4680.
      • 67. Sali A, Blundell TL (1993) Comparative protein modeling by satisfaction of spatial restraints. J Mol Biol 234: 779–815.
      • 68. Bastug T, Kuyucak S (2006) Energetics of ion permeation rejection binding and block in gramicidin a from free energy simulations. Biophys J 90: 3941–3950.
      • 69. Bastug T, Kuyucak S (2009) Importance of the peptide backbone description in modeling the selectivity filter in potassium channels. Biophys J 96: 4006–4012.
      • 70. Chen PC, Kuyucak S (2012) Developing a comparative docking protocol for the prediction of peptide selectivity proles: investigation of potassium channel toxins. Toxins 4: 110–138.
      • 71. Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, et al. (2005) Scalable molecular dynamics with NAMD. J Comput Chem 26: 1781–1802.
      • 72. Vanommeslaeghe K, Hatcher E, Acharya C, Kundu S, Zhong S, et al. (2010) CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force field. J Comput Chem 31: 671–690.
      • 73. Kumar S, Bouzida D, Swensen RH, Kollman PA, Rosenberg JM (1992) The weighted histogram analysis method for free-energy calculations on biomolecules. J Comput Chem 13: 1011–1021.
      • 74. Wilson MJ, Yoshikami D, Azam L, Gajewiak J, Olivera BM, et al. (2011) µ-Conotoxins that differentially block sodium channels NaV1.1 through NaV1.8 identify those responsible for action potentials in sciatic nerve. Proc Natl Acad Sci USA 108: 10302–10307.
      • 75. Chen R, Robinson A, Chung SH (2014) Mechanism of µ-conotoxin PIIIA binding to the voltage-gated N+ channel NaV1.4. PLoS One 9: e93267.