Voltage-gated sodium channels are targets for many drugs and toxins. However, the rational design of medically relevant channel modulators is hampered by the lack of x-ray structures of eukaryotic channels. Here, we used a homology model based on the x-ray structure of the NavAb prokaryotic sodium channel together with published experimental data to analyze interactions of the μ-conotoxins GIIIA, PIIIA, and KIIIA with the Nav1.4 eukaryotic channel. Using Monte Carlo energy minimizations and published experimentally defined pairwise contacts as distance constraints, we developed a model in which specific contacts between GIIIA and Nav1.4 were readily reproduced without deformation of the channel or toxin backbones. Computed energies of specific interactions between individual residues of GIIIA and the channel correlated with experimental estimates. The predicted complexes of PIIIA and KIIIA with Nav1.4 are consistent with a large body of experimental data. In particular, a model of Nav1.4 interactions with KIIIA and tetrodotoxin (TTX) indicated that TTX can pass between Nav1.4 and channel-bound KIIIA to reach its binding site at the selectivity filter. Our models also allowed us to explain experimental data that currently lack structural interpretations. For instance, consistent with the incomplete block observed with KIIIA and some GIIIA and PIIIA mutants, our computations predict an uninterrupted pathway for sodium ions between the extracellular space and the selectivity filter if at least one of the four outer carboxylates is not bound to the toxin. We found a good correlation between computational and experimental data on complete and incomplete channel block by native and mutant toxins. Thus, our study suggests similar folding of the outer pore region in eukaryotic and prokaryotic sodium channels.
Sodium channels, which play key roles in physiology of excitable cells, have been a subject of experimental studies during the last decades (Trimmer et al., 1989; Priestley, 2004; Wood et al., 2004; Rogers et al., 2006). The pore-forming α subunit of voltage-gated sodium channels (Nav1) is composed of four homologous repeats (I–IV), which are quasi-symmetrically arranged around the central pore. Each repeat contains six transmembrane segments (S1–S6), a membrane-diving P-loop (P), and the intracellular N and C termini. P-loops and extracellular thirds of transmembrane S5 and S6 helices contribute to the outer pore region, whereas cytoplasmic two-thirds of helices S5 and S6 contribute to the inner pore region.
Actions of naturally occurring toxins and synthetics drugs on Nav1 channels have been intensively studied. These include the hallmark sodium channel blockers tetrodotoxin (TTX), saxitoxin, and μ-conotoxins; sodium channel agonists batrachotoxin, veratridine, aconitine, and grayanotoxin; local anesthetics; and many other ligands (French et al., 1984; Moczydlowski et al., 1984; Linford et al., 1998; French and Terlau, 2004; Eijkelkamp et al., 2012; Waszkielewicz et al., 2013). However, in the absence of x-ray structures of eukaryotic sodium channels, understanding of atomistic mechanisms of action is still incomplete for some classes of ligands. X-ray structures of potassium channels in the closed and open states (Doyle et al., 1998; Li et al., 1998; Jiang et al., 2003; Long et al., 2005) have been used to build homology models of sodium channels, which provide structural explanation for numerous experimental data (Lipkind and Fozzard, 2000; Tikhonov and Zhorov, 2005; Choudhary et al., 2007; Tikhonov et al., 2014). However, potassium channels are poor templates for homology modeling of the outer pore region, which includes the selectivity filter. In potassium channels, the TVGYG motif downstream from the P-loop turn forms a narrow tunnel where five rings of backbone carbonyls catalyze the selective ion permeation. In sodium channels, the selectivity is provided by side chains of the D, E, K, and A residues from four P-loops (the DEKA ring), and the selectivity-filter region is wide enough to accommodate blockers as large as TTX. Therefore, attempts to rationalize the data on ligand action in the outer pore necessarily involved speculations about folding of the outer pore region (Lipkind and Fozzard, 2000; Tikhonov and Zhorov, 2005; Choudhary et al., 2007).
The x-ray structures of NavAb and other prokaryotic sodium channels have provided a new and possibly more reliable structural basis for understanding ligand action in the outer pore region (Payandeh et al., 2011; McCusker et al., 2012; Zhang et al., 2012). Selectivity in these channels is provided by a single EEEE ring. Unlike potassium channels, each subunit of a prokaryotic sodium channel has two helices before (P1) and after (P2) the selectivity filter glutamate. The P2 helices line the outer pore and form relatively wide extracellular vestibule. Action of TTX has been recently modeled in view of the NavAb folding (Tikhonov and Zhorov, 2012). However, TTX interacts only with a local region between the DEKA ring and the more extracellular ring of outer carboxylates. It remains unclear whether or not available structures of prokaryotic sodium channels can be used to rationalize action of much bigger toxins that target the outer pore.
Particularly interesting is the action of μ-conotoxins (Olivera and Cruz, 2001; Lewis et al., 2012). Conotoxins are biologically active peptides found in the venom of cone snails. μ-Conotoxins are rather small peptides (16–26 amino acids) with several disulfide cross-links that stabilize a typical folding pattern called framework III (Lewis et al., 2012). μ-Conotoxins do not target precisely the same binding area as TTX and saxitoxin (Zhang et al., 2010) and do not strongly interact with the selectivity filter, but interact with outer carboxylates and several residues downstream. The family of μ-conotoxins includes three toxins, which are particularly well studied: GIIIA (Cruz et al., 1989), PIIIA (Shon et al., 1998), and KIIIA (Zhang et al., 2007). GIIIA, the best-studied μ-conotoxin, binds with a high affinity to the Nav1.4 sodium channel (Sato et al., 1991; Becker et al., 1992; Dudley et al., 1995, 2000; French et al., 1996; Chang et al., 1998; Li et al., 2001, 2002; Cummins et al., 2002; Hui et al., 2002, 2003; Xue et al., 2003; Choudhary et al., 2007). PIIIA is a close homologue of GIIIA, but with some distinguishing features, which make it an interesting object for studying (Safo et al., 2000; McArthur et al., 2011a,c; Chen and Chung, 2012; Tietze et al., 2012). KIIIA, a short peptide of 16 amino acids (Zhang et al., 2007), is of special interest because of its small size and ability to bind to the channel simultaneously with TTX (Zhang et al., 2009; Van Der Haegen et al., 2011; Wilson et al., 2011; Khoo et al., 2012; Stevens et al., 2012).
Previous models of μ-conotoxin binding (Choudhary et al., 2007; McArthur et al., 2011b,c) were elaborated to rationalize solid experimental data using KcsA-based homology models. The fact that the outer pore regions in KcsA and NavAb have different folding calls for reconsidering structural aspects of μ-conotoxin–channel interactions. Here we used our NavAb-based model of Nav1.4 (Tikhonov and Zhorov, 2012) to dock μ-conotoxins GIIIA, PIIIA, and KIIIA in the outer pore region. The predicted complexes of Nav1.4 with μ-conotoxins are in agreement with a large body of experimental data. Furthermore, they allowed us to rationalize previously unexplained experimental data on residual currents in μ-conotoxin–bound channels.
MATERIALS AND METHODS
All calculations were performed by using the ZMM program (ZMM Software, Inc.). Nonbonded interactions were calculated with the AMBER force field (Weiner et al., 1984, 1986) and a cutoff distance of 8 Å. Electrostatic interactions were calculated by using the distance- and environment-dependent dielectric function (Garden and Zhorov, 2010). No specific energy terms were used for cation-π interactions, which were accounted for with partial negative charges at the aromatic carbons (Bruhova et al., 2008). Bond lengths and bond angles were kept rigid during the calculations.
The model of rat Nav1.4 was taken from our previous study (Tikhonov and Zhorov, 2012). The sequence alignment is shown in Fig. 1. Initial conformations of toxins were taken from the Protein Data Bank (PDB IDs are 1TCJ for GIIIA, 1R9I for PIIIA, and 2LXG for KIIIA). Voltage-sensing domains were not modeled because they are far from the binding site of μ-conotoxins. Big extracellular loops were truncated in our model because no experimental data on their folding, as well as interaction with μ-conotoxins, is available.
The Monte Carlo minimization (MCM) method (Li and Scheraga, 1987) was used to optimize the models. All side chain torsions were randomly sampled during the MCM trajectories. Both side chains and backbones were flexible during energy minimizations. MCM of each model was performed until 6,000 consecutive energy minimizations did not decrease the energy of the apparent global minimum. To incorporate experimental data in our calculations we used a system of constraints. A distance constraint is a flat-bottom energy function that allows a particular atom–atom distance to deviate energy-free between lower and higher limits and imposes a parabolic energy penalty if the distance is beyond the limits. Two types of constraints were used. To maintain the template folding we used “pin” constraints. A pin allows an α carbon of an amino acid residue to deviate up to 1 Å from the respective template position without a penalty and imposes a parabolic energy function to penalize larger deviations. Pins were applied to all residues in the P1 and P2 helices, but not to the selectivity filter region (residues 176–181 in NavAb). Pins are not suitable for maintaining the folding of conotoxins because the latter were mobile in the docking procedure. Instead we applied atom–atom distance constraints to keep distances between α carbons of a toxin close to values obtained in the respective nuclear magnetic resonance study. To bias experimentally identified toxin–channel contacts, distance constraints between centers of charged groups in the contacting amino acids were used. For all constraints, the energy penalty was calculated using the force constant of 10 kcal/mol/Å2. To minimize possible influence of constraints on the docking results, each optimization was performed in two stages. After MCM with constraints (the first stage), the model was refined by the second MCM, in which all constraints except pins were removed (the second stage).
Online supplemental material
Figs. S1, S2, and S3 show details of GIIIA docking procedure. Tables S1, S2, and S3 show residue–residue interaction energies for Nav1.4 complexes with GIIIA, PIIIA, and KIIIA, respectively. Figs. S4, S7, and S8 show main contacts for these complexes. Fig. S5 shows interactions that impose orientation of GIIIA in the model. Fig. S6 and Video 1 demonstrate switching of contacts between GIIIA and the channel because of the flexibility of long side chains of charged residues. Fig. S9 shows distributions of favorable ion positions in the model with wild-type and mutant toxins.
Selection of the model
Available x-ray structures of prokaryotic sodium channels (NavAb, NavMs, and NavRh) are very similar in the outer pore region. Particularly the P1 and P2 helices are very close at the superimposed x-ray structures (Fig. S1). Therefore, selection of a particular template is not essential for the modeling. In contrast, the sequence alignment between prokaryotic and eukaryotic channels is critical, as only one position shift in the alignments results in ∼100° reorientation of a residue in the α helix. The sequence similarity between prokaryotic and eukaryotic sodium channels is rather small. For example, the sequence identity between P-loops of NavAb and repeat I of Nav1.4 is only ∼21%. Recently, we demonstrated that a straightforward sequence alignment (without insertions/deletions) does not allow us to build a homology model consistent with the experimental data on TTX binding (Tikhonov and Zhorov, 2012). We proposed an adjusted alignment with insertions/deletions in the selectivity filter region (Fig. 1). The adjusted alignment allowed us to build a 3D model in which TTX-sensing residues are oriented inside the pore, whereas long polar residues before and after the selectivity filter form interhelical H bonds that stabilize the structure. The orientations of the outer carboxylates in the models with the straightforward and adjusted sequence alignments are significantly different (Fig. 2 A).
We compared the abilities of the NavAb-based models of Nav1.4 with the straightforward and adjusted sequence alignments to reproduce available data on GIIIA binding. Particularly valuable are results of the mutant cycle analysis, which revealed pairwise interactions between the toxin and the channel residues (Table 1). Here we used experimental data (Choudhary et al., 2007) to define a set of distance constraints, which would impose H bonds or salt bridges between side chains of respective residues in the channel and toxins. In some cases, unambiguous definition of specific pairwise channel–toxin contacts was impossible. For example, according to experimental data, outer carboxylates of repeats III and IV of the channel are involved in contacts with three basic residues of the toxin (Choudhary et al., 2007). Therefore, we have chosen only the strongest pairwise contacts for these residues (Table 1).
We placed GIIIA above the channel manually in such a way that the toxin mass center was at the pore axis and the toxin face, which contains basic residues (K9, K11, R13, K16, and R19), was oriented toward the channel. The distance between the plane of α carbons of the outer carboxylate ring and the toxin backbone atom, which is closest to the ring, was 11 Å. At this distance, toxin did not interact with the channel. Using this general orientation, we generated 36 starting toxin–channel mutual dispositions by consecutively turning the toxin around the pore axis. From each starting point, the energy was Monte Carlo–minimized with the toxin–channel distance constrains. A representative MCM trajectory is shown in Fig. S2. Next, each model was refined by an unconstrained MCM. We obtained generally the same orientation of the toxin for every starting point in models with both the straightforward and adjusted alignments. Fig. S3 shows the superimposition of all the models obtained from different GIIIA orientations for the model with the adjusted alignment. This result indicates that the system of toxin–channel distance constraints unambiguously determines orientation of the toxin.
In both models, the large toxin body did not fit into the selectivity filter. Rather it bound above the narrowing vestibule. Side chains of GIIIA basic residues reached the outer carboxylates, which are located at the extracellular entrance to the vestibule. Importantly, the toxin mass center shifted from the pore axis toward repeat III (Fig. 2 E). The C-end helix of the toxin strongly interacted with the P2 helix of repeat III. The N-end of the toxin projected to the extracellular space. The middle part, which contains arginine R13, occurred within the outer pore, near repeat II.
Despite this general similarity, results of docking showed significant advantage of the model with the adjusted alignment. For this model, the root-mean-square deviation (RMSD) of α carbons in the toxin-binding region from the NavAb structure was 1.39 Å. During the refining unconstrained MCM, the toxin slightly moved (the RMSD of α carbons between the constrained and unconstrained models was as small as 0.5 Å). In contrast, the experimental distance constraints caused significant deformation of the model, which is based on the straightforward sequence alignment. Indeed, the RMSD for α carbons in the toxin-binding region between the model and the NavAb structure was as large as 2.04 Å. The refining unconstrained MCM shifted the toxin ∼4.5 Å away from the channel (Fig. 2 B). The toxin–channel interaction energy in the straightforward alignment model was 33% higher (less preferable) than in the adjusted alignment model.
We calculated energy of interactions between individual residues of the toxin and the channel in the two refined models (Fig. 2, C and D). The toxin–channel pairwise interactions were plotted against the experimental free energies of respective contacts (Table 1). Note that not all of the experimental contacts were used as constraints during the docking procedure. For the straightforward alignment model, the correlation is weak (Pearson correlation coefficient P = −0.44 and significance of correlation P = 0.17). Contrarily, for the adjusted alignment model, P = −0.71 and the correlation appeared statistically significant (P = 0.02).
These results clearly show that the adjusted alignment model fits experimental data much better than the straightforward alignment model. Therefore, all further analysis was performed with the adjusted alignment model.
Binding of GIIIA
Although our model readily reproduced many experimentally found specific toxin–channel interactions (Table S1 and Fig. S4), the contact between hydroxyproline hP17 and M1240 (Dudley et al., 2000) was not observed. hP17 leaned toward repeat III P2-helix (DIIIP2), where M1240 is located, but no specific interactions between the two residues were found. The hydroxyl group of hP17 occurred in a close proximity to D1248 at the C-end of DIIIP2. We imposed an additional distance constraint to bias an H bond between the hydroxyl group of hP17 and carboxyl oxygen of D1248 and obtained a model (Fig. S5, top) in which the contacts obtained in the previous model were preserved. Thus, our model provides a structural explanation for the role of hP17 in the toxin action. Importantly, this new contact, which is predicted by our calculations, is not yet reported in experimental studies. D1248 does not have any matching residues in homologous positions of other repeats, and thus the new toxin–channel contact may be an important determinant of the toxin–channel mutual disposition.
Substitutions of D762 and E765 in DIIP2 affect GIIIA binding (Xue et al., 2003). In our model these residues are located far from GIIIA and cannot directly interact with it. However, according to the NavAb crystal structure, D762 stabilizes the outer pore by making a salt bridge with R395. E765 is situated on the same side of the P2-helix as D762, thus facing toward the P1-helix of repeat I. According to our model, both D762 and E765 make specific contacts with R395, contributing to stabilization of the outer pore folding (Fig. S5, bottom). We suggest that alanine substitutions of D762 and E765 (Xue et al., 2003) affect the toxin binding allosterically rather than directly.
In our model the body of GIIIA, which contains the α helix, interacts mainly with the P2 helix of repeat III. Interestingly, repeat III aspartate D1248, which is highly conserved in sodium channels, is an attractive target for electrostatic interaction with the toxin. Other repeats of the channel have either neutral or positively charged residues in the matching positions of the sequence alignment (Fig. 1). Another possible cause for the asymmetric binding of the toxin is the asymmetric ring of the outer carboxylates. Although the outer carboxylates in repeats I, II, and IV are in the matching positions of the alignment, repeat III D1241 is shifted one position toward the C-end.
The calculated energies of pairwise interactions correlate reasonably well with the experimentally determined values. However, rather poor correlation was found for some medium-strength contacts between the outer carboxylates in repeats III and IV and toxin basic residues (Fig. 2 D). A possible reason is that we considered only a single-complex structure and thus our calculations ignored the entropy. Because of the large conformational flexibility of long side chains in the toxin basic residues and the channel acidic residues, the experimentally determined pairwise energies likely integrate contributions from many micro-states of the channel–toxin complex. To explore this possibility, we performed additional calculations. We used the above described model to generate 20,000 starting points with random conformations of side chains of GIIIA basic residues. From each starting points the energy was Monte Carlo–minimized with the backbone geometry preserved by the system of distance constrains between α carbons of the toxin (see Materials and methods). The results are shown in Fig. 2 F. In the ensemble of lowest-energy complexes, the charged terminal groups in flexible side chains of basic residues fluctuated up to 5 Å. This flexibility may switch some contacts (Fig. S6 and Video 1). Key basic residues of the toxin are K11, R13, K16, and R19. According to our results, K11 mostly interacted with the channel repeats I, III, and IV, but not with repeat II because of the presence of R13 in close proximity of repeat II. R13 interacted mostly with repeats I and II, but can interact with any repeat because of its long side chain. K16 in the middle of the helical part of the toxin can interact with any channel repeat, depending on contacts of other toxin residues. R19 can only reach repeats III and IV because it is situated close to the C-end of the toxin, downstream from the helical part, and cannot approach closer repeats I and II. Thus, our model predicts that within the same binding mode and backbone geometry of the toxin, its basic residues can form different contacts with acidic residues of the channel.
Binding of PIIIA
Residues R12, R14, K17, and R20 are important for PIIIA–channel interactions (McArthur et al., 2011c). However, pairwise contacts between PIIIA and Nav1.4 are not known. Therefore, the docking methodology, which we used for GIIIA, is not applicable for PIIIA. Because PIIIA folds like GIIIA, we suggested that in the complex with the channel, the α helix and key basic residues of PIIIA are oriented like those of GIIIA and generated respective starting point for our calculations. The Monte Carlo–minimized binding mode is shown in Fig. 3. In our model, the key PIIIA residues are involved in the following pairwise interactions: R12/E403, R14/E758, K17/D1241, K17/D1532, and R20/D1532 (Fig. S7 and Table S2).
Although the overall binding modes for GIIIA and PIIIA are rather similar (RMSD for α carbons of key basic residues is 1.8 Å), there are certain differences in pairwise toxin–channel interactions. For example, specific interaction between D12 and T759 (Choudhary et al., 2007) significantly contributed to the stability of the GIIIA–channel complex in our model. PIIIA lacks an acidic residue in the homologous position. As a result, in our models, specific interactions of arginine R12 in PIIIA differ from interactions of the homologous lysine K11 in GIIIA. R12 in PIIIA interacted strongly with repeat I, whereas K11 in GIIIA interacted with repeats III and IV because contacts of K11 with repeat I were sterically precluded by D12. R14 formed tight contacts with repeats I and II, but weakly interacted with repeat III, and even weaker with repeat IV. K17 interacted with repeat III and IV like its analogue in GIIIA, K16. Finally, R20 (the analogue of R19 in GIIIA) interacted only with repeat IV. Like in the case of GIIIA–channel interactions, the shown macrostate may represent many microstates that exist as a result of the possibility of switching contacts between long flexible side chains.
McArthur et al. (2011c) explored the voltage dependence of action of PIIIA and its mutants and determined the depths (δ) of residues in the membrane electric field. The δ values rank as follows: R14 > K17 > R20 > R12 > S13 > R2 > G6. In our model, the relative depth of the charged residues (determined by positions of their basic atoms) is ranked as follows: R14 > K17 > R20 > R12 > R2. S13 is located approximately at the same depth as R12, whereas G6 and R2 are located at the extracellularly projected N-end of the toxin. Thus, the toxin-binding orientation in our model is consistent with experimental data. It should be noted that no data on voltage dependence were used to obtain the toxin orientation.
Binding of KIIIA
KIIIA, the smallest μ-conotoxin, lacks an analogue of the N-terminal segment, which is present in GIIIA and PIIIA, and contains only three basic residues in the α helix (Fig. 1). The data on pairwise interactions (McArthur et al., 2011b) suggest that all of these residues interact with repeat III outer carboxylate (D1241). We used these data to impose distance constraints, which clamped H bonds between the toxin basic side chains and the outer carboxylate of repeat III. MCM yielded a model (Fig. 4) in which the position and orientation of the KIIIA α helix are similar to those of PIIIA and GIIIA. The entire set of contacts of basic residues in KIIIA is shown in Fig. S8 and Table S3. Importantly, because of its small size and the lack of N-terminal segment, which would be analogous to N-terminal segments in GIIIA and PIIIA, KIIIA did not cover the pore with its body and weakly interacted with repeat I.
Intriguingly, TTX can reach its binding site deeply in the outer pore in the presence of KIIIA, although with a slower kinetics than in the absence of KIIIA (Zhang et al., 2009). To rationalize these data, we have built a ternary complex model of the channel with KIIIA and TTX. We used the TTX-binding mode from Tikhonov and Zhorov (2012) and KIIIA-binding mode described above. The ternary complex model (Fig. 4, A and B) did not reveal toxin–channel or toxin–toxin sterical clashes. The compact TTX molecule is bound deeply in the narrow portion of the pore between the DEKA ring and the ring of outer carboxylates, whereas the larger KIIIA toxin did not penetrate deep into the pore, approached the outer carboxylates from the extracellular side, and formed salt bridges with them. Each carboxylate can accept several H bonds so that simultaneous contacts with TTX and KIIIA were easily formed. The repeat III outer carboxylate, which unlike other outer carboxylates faces away from the pore axis, did not interact with TTX, but formed a salt bridge with R10 of KIIIA.
To find a possible passage for TTX between KIIIA and the channel, we pulled TTX out from the ternary complex model in the extracellular direction. To do this we imposed a distance constraint between the central carbon atom in the TTX guanidinium group and a plane drawn through four α carbons of the DEKA ring. Next we systematically increased the minimal value of the distance constraint with the step 1 Å and Monte Carlo–minimized energy in each step. The lowest-energy structure found in a given step was used as the starting point for the next step. During the calculations, the binding mode of KIIIA was preserved by distance constraints, which clamped contacts of its basic residues with the outer carboxylates. The TTX egress trajectory is rather smooth (Fig. 4 C). The corresponding energy plot (Fig. 4 D) does not have energy barriers, indicating the absence of steric obstacles in the predicted pathway for the TTX egress.
Mechanisms of incomplete block
Experimental studies revealed incomplete block of Nav1 channels by KIIIA and some GIIIA and PIIIA mutants (Hui et al., 2002; McArthur et al., 2011a; Wilson et al., 2011), but structural rationale for these data is lacking. To address the problem, we used the approach that was used to predict the TTX egress pathway. A sodium ion was placed in the DEKA ring, and the egress pathway in the presence of GIIIA was calculated. The ion trajectory is interrupted near the toxin (Fig. 5 A). The energy profile contains zones of negative values, where the ion experiences favorable electrostatic interactions with the DEKA ring and the outer carboxylates, and the zone of positive values, which corresponds to repulsive interactions with the toxin (Fig. 5 C).
Next we computed a similar energy profile for the R13A mutant of GIIIA bound in Nav1.4. We did not search systematically the mutant-binding mode, as we did for GIIIA, but just replaced R13 with alanine in our GIIIA–channel model and then Monte Carlo–minimized the energy of the complex. The toxin body did not move, and all of the toxin–channel contacts (except those at the mutated reside) were preserved. However, the egress pathway of the sodium ion in this model was quite different from that predicted for the wild-type toxin in terms of energy and geometry (Fig. 5 B). The trajectory became smooth. Comparison of panels A and B in Fig. 5 demonstrates a key role of mutation R13A in this difference. The zone of positive energies disappeared (Fig. 5 C), indicating that the mutant toxin does not impose obstacles for the ion movement between the selectivity filter and the extracellular space. Thus, our calculations allow us to explain the different action of GIIIA, which produces complete block, and its R13A mutant, which produces only a partial block.
To further explore possibilities of sodium movement through the toxin-bound channels, we generated 10,000 staring structures in which the ion was randomly placed in the outer pore region and briefly optimized the starting positions (10 minimizations in each of the 10,000 MCM trajectories). For the analysis, we selected only those structures in which the energy of interaction of the ion with the rest of the system was negative (attracting interactions dominate). This approach allowed us to find regions in the outer pore that can be populated by ions. The calculations were performed for the toxin-free channel model and for the complexes of the channel with the wild-type and mutant toxins. The results are shown in Fig. S9, and some representative data are shown in Fig. 6.
In the toxin-free Nav1.4 model, continuous occupancy of the outer pore region is seen (Fig. 6 A). A large number of energetically possible ion positions in the extracellular half of the outer pore reflects the large free space at this region. The bound GIIIA toxin displaces most of ions from this region (Fig. 6 B), but two constellations of ions are seen at the extracellular part of the outer pore and at the selectivity filter. Importantly, between these constellations there is a zone, which practically lacks energetically favorable positions for the ion (Fig. 6, B and E). It should be noted that the absence of favorable ion positions in this zone is not caused by complete sterical occlusion of the pore. Rather, the ions cannot occupy this region because of electrostatic repulsion with positively charged residues of GIIIA. Indeed, when we disabled the electrostatic interactions for the randomly placed ion, energetically favorable ion positions appeared all along the pore axis without a break (not depicted). Thus, in our model the nature of the channel block by μ-conotoxin is predominantly electrostatic.
For the R13A mutant of GIIIA, the shape of ion distribution is more similar to that in the toxin-free channel. Although many possible ion positions in the outer pore of toxin-free channel are eliminated by the toxin, such positions are seen at all the levels of the outer pore, implying that the pathway for the ions is not interrupted (Fig. 6 C). Similar uninterrupted ion distributions were obtained for KIIIA and PIIIA mutants, which produce only incomplete block. In contrast, for the GIIIA mutants (K9A, K11A, K16A, and R19A) and for the R20A mutant of PIIIA, which produce complete block (McArthur et al., 2011a), the ion distributions were similar to that of GIIIA wild type, i.e., two regions populated by ions are separated by the zone, where ions are repelled by the toxin (Fig. S9).
To analyze the obtained results, we generated a list of specific contacts in the channel complexes with the wild-type and mutated toxins (Table 2). We found that uninterrupted ion distributions in the models are observed if and only if at least one of the outer carboxylates is not involved in the interactions with the toxin basic residues. Thus, KIIIA, which has only three basic residues, is unable to block permeation completely. GIIIA and PIIIA, which completely block the current, specifically interact with all four outer carboxylates. The models with toxin mutants, in which any of the outer carboxylates is not engaged in specific interactions with the toxin, correspond to the incomplete block. If an outer carboxylate can form salt bridges with two or more basic residues of the toxin, mutation of a single basic residue does not release the outer carboxylate, which remains neutralized and therefore unattractive for a sodium ion. The block remains complete in such cases. For example, alanine substitution of R19 in GIIIA does not produce a residual current (McArthur et al., 2011a) because K16 and K11 readily interact with D1532, which is the strongest contact for the native R19 in GIIIA.
Thus, our study proposes a critical role of the outer carboxylates in the process of ion permeation. These residues likely serve as intermediate binding sites for ions moving from the extracellular space toward the selectivity filter (Fig. 6 D). If at least one of the outer carboxylates remains free, the ion permeation is possible, although it is reduced to some extent. This model agrees with experimental data (McArthur et al., 2011a,b) and with the previous modeling suggestion about the role of outer carboxylates (Tikhonov and Zhorov, 2007).
Effects of R13 substitutions in GIIIA were systematically studied (Hui et al., 2002). It was demonstrated that the value of residual current depends on both size and charge/polarity of a substitute. To reproduce these results, we modeled R13 substitutions by K, Q, W, Q, N, A, E, and D and calculated the ion distributions as described above. Ion distributions for all the mutants are not completely interrupted, suggesting the existence of some residual current (Fig. 6 E). To obtain quantitative correlation with experimental data (Hui et al., 2002), we plotted the number of ion at axis coordinate z = 8 Å (where the ion distribution is interrupted for the wild-type toxin) against the experimental data on residual current. Results of our calculations (Fig. 6 F) demonstrate a good correlation with experiments (R = 0.96, P = 3.4 × 10−5). Thus, the model reproduced residual currents in the channel complexed with different GIIIA mutants not just qualitatively, but quantitatively.
In the present work we have used a NavAb-based homology model of the outer pore region of the voltage-gated sodium channel Nav1.4 to dock μ-conotoxins GIIIA, PIIIA, and KIIIA and some of their mutants. A key question addressed in our work was whether or not the NavAb structure is suitable for rationalizing a large body of available experimental data on the action μ-conotoxins in eukaryotic channels. Unlike previous models of the outer pore region in Nav1 channels, which were necessarily based on potassium channel structures, we did not modify the NavAb template folding to obtain agreement between the modeling predictions and experimental data. Our calculations have demonstrated that the NavAb-based model does allow us to rationalize various experimental data. This result indicates that the outer pore in prokaryotic and eukaryotic voltage-gated sodium channels folds similarly. In other words, our results demonstrate that the NavAb structure allows us to proceed from “hand-made” models of the outer pore region in voltage-gated sodium channels to structure-based models.
This conclusion is not a trivial one. The sequences of NavAb and eukaryotic channels in the selectivity filter and outer pore regions are significantly different. NavAb contains the EEEE ring in the selectivity filter, whereas eukaryotic channels have the DEKA ring (Fig. 1). The EEEE ring is typical for calcium-selective rather than sodium-selective channels. Furthermore, NavAb lacks the ring of outer carboxylates, which has been demonstrated to participate in the ion permeation. This suggests that mechanisms of selectivity and permeation in NavAb and eukaryotic sodium channels have different structural bases.
Recently, PIIIA was predicted to block NavAb channel in the binding mode that significantly differs from that in our model of PIIIA with Nav1.4 (Chen and Chung, 2012). The large difference between sequences of prokaryotic and eukaryotic sodium channels in the outer pore region may explain different binding modes of μ-conotoxins in these channels. A model of PIIIA in Nav1.4 was proposed, which is based on the straightforward alignment of Nav1.4 and NavAb sequences (Chen et al., 2014). Using MD simulations and constraints that pulled the toxin basic side chains toward the selectivity filter, the authors obtained models in which PIIIA is located in the outer pore deeper than in our model. The proposed toxin–channel interactions also differ significantly from those in our models. In the absence of experimental data on PIIIA–channel pairwise contacts, it is hardly possible to select the most realistic PIIIA-binding model. It should be noted that we addressed diverse experimental data on the action of GIIIA, PIIIA, and KIIIA and the obtained models readily reproduced various experimental observations, which were not used during docking procedure. For instance, the binding mode elaborated for GIIIA explains simultaneous binding of KIIIA and TTX. The PIIIA-binding mode agrees with the experimental data on relative electric depths of individual residues.
Incomplete block by some μ-conotoxins and their mutants is well known (McArthur et al., 2011a,b), but structural bases of this effect were unclear. It is demonstrated that the current decreases because of decrease of the channel conductance (Wang et al., 2000, 2001). Our modeling suggests a molecular mechanism of this effect. Unlike TTX and related small-size ligands, μ-conotoxins do not penetrate into the narrow portion of the outer pore and do not completely occlude the pore lumen by their bodies. According to our models, conotoxins block the current by catching and neutralizing the outer carboxylates. If all four outer carboxylates are involved in the interactions with a toxin, the permeation is completely blocked. If some of the outer carboxylates remain free, the toxin only partially blocks the current. This proposed scheme (Fig. 6 D) agrees with the data that mutations of the outer carboxylates reduce the channel conductance. Although our models emphasize the role of electrostatic interactions in compete or incomplete channel block, the value of residual current obviously depends on the steric factors. This is particularly true for R13 substitutions in GIIIA because this residue directly faces the pore. The size of the side chain in this position largely determines the number of ions, which are displaced from the critical region of the outer carboxylates, in agreement with experimental data (Hui et al., 2002).
Based on our models, some testable predictions can be proposed. Particularly, we found that the hydroxyl group of hP17 can specifically interact with carboxyl group of D1248. Also, according to our model, D762 and E765 make specific contacts with R395 and contribute to stabilization of the outer pore folding, suggesting that alanine substitutions of D762 and E765 (Xue et al., 2003) affect the toxin binding allosterically rather than directly. From this point of view, other pairs of specifically interacting residues (Q407-E1524, Y1244-R750, N1536-Q1232), which stabilize mutual disposition of P1 and P2 helices, can also have a certain impact on conotoxin binding.
Certain limitations of our models should be discussed. First, we modeled the asymmetric channel using a symmetric template. Certainly, sequence differences between individual repeats should cause some structural asymmetry on the backbone. To take this effect into account, we did not clamp the template backbone but used pin constraints, which allow penalty-free deviation of α carbons up to 1 Å from the respective template positions. However, the asymmetry of the pore domain may be more significant, in particular because of the influence of the channel regions, which were not considered in our model. Second, we did not take into account all possible factors that affect conotoxin binding. For example, the toxin–channel affinity significantly depends on the ionic strength (Li et al., 2003). Qualitatively, it can be explained by competition between metal ions and basic residues of the toxin for interactions with acidic residues of the channel. However, precise quantification of this effect is difficult. Third, we did not perform hands-free docking of the toxins by systematically exploring all possible binding modes. Inherent limitations of the homology modeling, some of which are mentioned above, make hands-free calculations not completely reliable. That is why we used available experimental data as constraints for the docking procedure. Fourth, we assumed that point mutations do not change the overall binding mode. For some mutants this is not evident. For example, the R13A mutant of GIIIA subtly, but clearly stands apart from other mutants (Hui et al., 2002), suggesting that the binding mode is changed. However, at the present level of model precision, it is hardly possible to correctly predict the difference between the binding modes of wild-type and mutant toxins.
Given these limitations, we did not attempt to predict binding affinities of the native and mutant toxins. However, the limitations do not affect the main conclusions of our study. Our results suggest that NavAb-based models of the outer pore region of eukaryotic voltage-gated sodium channels provide significant progress in understanding the structural organization of these important proteins. Until the x-ray structure of a eukaryotic sodium channel becomes available, the NavAb-based model can be used to rationalize experimental data and design new experiments.
Computations were made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET). We thank an anonymous reviewer of the manuscript for bringing our attention to results of systematic analysis of GIIIA_R13 mutations and to the problem of steric versus electrostatic nature of the sodium channel block by μ-conotoxins.
This work was supported by grants to D.B. Tikhonov (the program of Russian Academy of Sciences “Molecular and Cellular Biology” and RFBR-13-04-00724) and to B.S. Zhorov (GRPIN-2014-04894 from the Natural Sciences and Engineering Research Council of Canada and RFBR-14-04-378).
The authors declare no competing financial interests.
Kenton J. Swartz served as editor.