INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Increasing efforts are devoted to the structural and functional characterization of K⁺ channels (Blaustein et al., 2000; Cha and Bezanilla, 1998; Doyle et al., 1998; Hong and Miller, 2000; Monks et al., 1999; Sokolova et al., 2001; Tiwari-Woodruff et al., 2000; Zhou et al., 2001). In this regard, the toxins isolated from various venoms, which are potent inhibitors of K⁺ channels, are of particular interest. These toxins are small proteins that block the passage of K⁺ ions by binding at the pore entryway on the extracellular side of the channel, thereby inhibiting the ion flux. The interactions with K⁺ channels are among the strongest and most specific known in protein-protein complexes. In virtue of their remarkable properties, these toxins can serve as powerful investigative tools in experimental studies of ion channels.

Charybdotoxin (CTX) and its close homolog Agitoxin (AgTx2) isolated from the venom of scorpion Leirus quinquestriatus are the best characterized channel blocking toxins. Their three-dimensional (3D) structures in solution have been determined to high resolution using nuclear magnetic resonance (NMR) (Bontems et al., 1991, 1992; Krezel et al., 1995). With dissociation constants in the nanomolar range, CTX is a potent inhibitor of the high conductance Ca²⁺-activated channels (Anderson et al., 1988) and of the voltage-activated Shaker K⁺ channels (Stocker and Miller, 1994), whereas AgTx2 is a strong inhibitor of Shaker (Miller, 1995). Other well-characterized channel blocking toxins are, to name a few, the voltage-gated K⁺ channel inhibitors ShK (Tudor et al., 1996) and BgK (Dauplais et al., 1997) extracted from sea anemones; tertiapin Q, a derivative of honey bee toxin blocker of the inward rectifier (Kir1) ROMK (Jin et al., 1999); and dendrotoxins isolated from mamba Dendroaspis snake venoms (Harvey, 1997) such as -dendrotoxins, a potent blocker of Kv1.1 and Kv1.2 (Hurst et al., 1991; Stocker et al., 1991; Tytgat et al., 1995), and -dendrotoxins, a blocker of the inward rectifier Kir1 (Imredy et al., 1998) as well as voltage-activated channels (Dolly and Parcej, 1996; Pongs, 1990; Pongs, 1993).

The scorpion toxins CTX and AgTx2 have been used to identify the pore region of K⁺ channels (MacKinnon et al., 1990; MacKinnon and Miller, 1989; Rauer et al., 2000), determine the channel subunit stochiometry (MacKinnon, 1991), elucidate the structural topology of the extracellular surface of channels (Goldstein et al., 1994; Hidalgo and MacKinnon, 1995; Ranganathan et al., 1996), and assign the functional sidedness of the proton activation of KcsA (Heginbotham et al., 1999). Their interaction with K⁺ channels has been the object of numerous studies (Gross and MacKinnon, 1996; Hidalgo and MacKinnon, 1995; MacKinnon and Miller, 1988; Miller, 1995; Naini and Miller, 1996; Naranjo and Miller, 1996; Ranganathan et al., 1996; Rauer et al., 2000).

Using the known 3D structure of the toxins as a template (Bontems et al., 1991, 1992; Krezel et al., 1995), thermodynamic analysis of double mutant cycles with CTX (Naini and Miller, 1996; Naranjo and Miller, 1996; Rauer et al., 2000; Stocker and Miller, 1994) and AgTx2 (Gross and MacKinnon, 1996; Hidalgo and MacKinnon, 1995; Ranganathan et al., 1996) was used to extract spatial restraints between pairs of residues (one on the toxin and one on the channel) and structurally characterize the pore of Shaker. As shown in the case of the Barnase-Barstar complex, for which a crystallographic structure of the complex is available (Buckle et al., 1994; Guillet et al., 1993), strong couplings in double mutant cycles are, on average, generally correlated with short distances between the mutated residues in the protein-protein complex (Schreiber and Fersht, 1995). The subsequent determination of the structure of the KcsA K⁺ channel using x-ray crystallography (Doyle et al., 1998) made it possible to verify and confirm the structural information deduced from the analysis of double mutant cycle (MacKinnon et al., 1998). Nonetheless, although this initial analysis established the structural similarity of KcsA to Shaker and provided many of the essential elements for understanding the microscopic basis for the specificity of toxin-channel interactions, several aspects require additional considerations. In particular, only the largest toxin-channel contacts were included in the analysis and a large amount of the information available was not used. In addition, it was implicitly assumed that the coupling constants could be translated into interresidue distances, which is true only on average. Thus, the structure of the toxin-channel complex has not been determined at the highest possible resolution that can be achieved with the currently available information from the double mutant data (Hidalgo and MacKinnon, 1995; Ranganathan et al., 1996).

To make progress in understanding the microscopic basis for the specificity of these toxins, it is necessary to characterize the interactions in the toxin-channel complex at the atomic level. There is, however, currently no high-resolution 3D structure available from experiments of any toxins in complex with a K⁺ channel. The goal of this study is to determine the 3D structure of AgTx2 in complex with the extracellular vestibule of the Shaker K⁺ channel at atomic resolution using computational methods and available experimental data. A particular emphasis is put on assessing the quality of the models with some degree of confidence. The computational protocol can be divided into two parts. First, a large number of protein-protein complexes are generated using the distance restraints deduced from thermodynamic double mutant data. Second, the quality of the generated models is examined by generating molecular dynamics (MD) trajectories of the toxin-channel complex solvated by explicit water molecules and by calculating the relative binding free energy of various mutant complexes using a continuum solvent model and comparing with available experimental values. The docking protocol, allowing full flexibility of the side chains, differs from the rigid-body Brownian dynamics approach of Cui et al. (2001) who have simulated the diffusional encounter of Lq2 with KcsA. The present MD treatment is more similar to that used for docking CTX and ShK to different voltage-gated K⁺ channels (Kalman et al., 1998; Rauer et al., 2000).

The work is presented as a computational method aimed at extracting structural information from experimental data. In this regard, the outcome differs from pure modelization of protein-protein complexes (Camacho and Vajda, 2001; Palma et al., 2000; Ritchie and Kemp, 2000). In the next section, the computational methods that are used for the structure determination and validation are described. The methods are first tested on the Barnase-Barstar complex, for which both a high-resolution x-ray structure (Buckle et al., 1994; Guillet et al., 1993) as well as plenty of experimental data from thermodynamic double mutant cycles are available (Schreiber and Fersht, 1995). Thereafter, we describe the structure determination of the AgTx2/Shaker complex as well as the evaluation of the resulting binding modes and the results are discussed. Finally, the paper is concluded with a summary of the most important findings.

	THEORY AND METHODS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Modeling of Shaker

Although there is a growing body of information about the structure of Shaker, its 3D structure at atomic resolution remains unknown. At the present time, the only ion channel for which a 3D structure at atomic resolution is available are the KcsA channel from Streptomyces lividans (Doyle et al., 1998; Zhou et al., 2001), and the MthK channel from Methanobacterium thermoautotrophicum (Jiang et al., 2002). Both Shaker and KcsA are made of four identical subunits disposed symmetrically around a common axis corresponding to the pore. Analysis of the amino acid sequence suggests that each subunit of Shaker contains six putative transmembrane (TM) segments called S1 to S6 (Jan and Jan, 1997; Tempel et al., 1987). The S1 to S4 segments are involved in the voltage-gating mechanism (Jan and Jan, 1997; Liman et al., 1991; Logothetis et al., 1992), whereas the S5 and S6 segments form the central pore region responsible for the conduction and selectivity of K⁺ ions (Heginbotham et al., 1992, 1994; Zhou et al., 2001). Although the subunits of KcsA comprises only two transmembrane helices (TM1 and TM2), the amino acid sequence is very similar to the segment S5-S6 forming the central pore of Shaker: 31 of 93 residues are identical and 15 are similar for a global sequence similarity of ~49% (Cortes and Perozo, 1997; Doyle et al., 1998; Schrempf et al., 1995). In addition, the structural similarity of Shaker relative to KcsA is also supported by experiments: AgTx2 binds to a triply mutated version of KcsA (although more weakly than to Shaker), indicating that the two channels do share the same pore structure (MacKinnon et al., 1998). Furthermore, a chimera incorporating the pore domain of KcsA into Shaker was shown to retain the main functional features of Shaker (Lu et al., 2001). All these observations suggest that a reliable comparative model of the central pore of Shaker can be constructed using the crystallographic structure of the KcsA channel as a template (Fiser et al., 2000; Sali et al., 1995).

The similarity-based models of Shaker were generated by optimizing a larger number of internal spatial restraints extracted from the KcsA structure using the program MODELLER v. 4 (Sali and Blundell, 1993). This technique has been shown to yield 3D models that are generally accurate (Sali and Blundell, 1993). One-hundred 3D models of Shaker were generated. The ideal fourfold symmetry of the tetrameric channel was imposed on all the models. The best model (the one with the lowest MODELLER restraint energy), shown in Fig. 1 b, was kept and then further refined with energy minimization of the side chains. Two potassium ions were inserted at positions corresponding to the upper inner site (in contact with the carbonyl of Thr-75 and Val-76 in KcsA) and the cavity site of KcsA (Doyle et al., 1998). The outermost sites were left unoccupied as it can be anticipated that they will clash with the Lys-27 side chain of the toxin (see below). The segment S6 was modeled on the basis of the crystallographic structure despite indications that its structure may be different toward the intracellular end beyond residue Val-474 (Holmgren et al., 1998). However, these structural differences in S6 are of no consequence on the binding of AgTx2 to the extracellular side of the central pore and can be safely ignored in the current modeling.

View larger version (34K): [in this window] [in a new window]   FIGURE 1   (a) Chosen NMR structure of AgTx2 (Krezel et al., 1995). Residues Arg-24 and Lys-27, which are important for the binding to Shaker, are colored blue. (b) Comparative model of the Shaker K+ channel (S5-S6). The front and rear monomers have been removed for easier visualization. Coloring corresponds to the following regions: dark blue, outer helix (S5, 391-417); gray, turret (418-424); orange, pore helix (425-440); blue, inner helix (S6, 454-482). The backbone of the selectivity filter is shown explicitly. Monomers are labeled A, B, C, and D in the clockwise sequence from the extracellular side.

Docking and simulated annealing

Coupling constants () derived from thermodynamic analysis of double mutant cycles represent the main source of information about the structure of the toxin-channel complex. It has been shown that the magnitude of  reflects the spatial closeness of two residues in the protein-protein complex (Schreiber and Fersht, 1995). The pairs of residues exhibiting the strongest coupling in the thermodynamic analysis of double mutant cycles for AgTx2 and Shaker are given in Table 1 (Ranganathan et al., 1996). Typically, large couplings are indicative of a short distance between the mutated residues. In contrast, weak couplings are slightly less informative because they may reflect a large distance but could also arise from a cancellation of effects. Although such weak couplings can still be useful in postanalysis to validate the protein-protein complex, it is therefore preferable to use only the strongest couplings as distances restraints in generating a 3D model of the protein-protein complex.

In the case of the AgTx2-Shaker complex, one practical problem arises from the fact that the experimentally observed couplings correspond to an average over four possible and equivalent orientations of the toxin bound to the tetrameric channel. For this reason, a distance restraint cannot be assigned unambiguously between a residue in the toxin and a residue in a specific channel subunit. To resolve the intrinsic ambiguity from the fourfold symmetry of the system, the distance restraint between a residue of the toxin (i) and of the channel (j) were implemented as

(1)

in which r is the upper bound on the distance restraint, K is the strength of the restraint and
In the case of polar or charged residues, the minimum over all distances between the heavy atoms of the polar or charged groups was taken, whereas for nonpolar residue pairs the minimum over all distances between all heavy atoms of each side chain was taken. Such spatial restraint is "ambiguous" because it is not assigned between a specific pair of atoms (one in the toxin and one in a given subunit of the channel) but is effectively applied only between the toxin residue and the residue of the channel subunit (A, B, C, or D) that is the closest instantaneously. Standard distance restraints between pairs of specific residues were applied in the later stages of the structural refinement as the contacts between the toxin and the channel were unambiguously characterized.

The toxin was docked onto the channel surface in the presence of a selected set of distance restraints taken from Table 1. The docking protocol consists of a series of short MD trajectories followed by energy minimizations. The channel was positioned such that the wide nonpolar cavity was centered at Z = 0, which corresponds roughly to the center of a membrane bilayer, and the pore was aligned with the Z-axis with the extracellular side on the positive side. The atoms of all the residues of the channel for which Z < 9 Å (corresponding to residues below Gly-444 in the selectivity filter) were kept fixed in all the MD calculations. For the remaining residues of Shaker, the backbone atoms were restrained relative to the initial model using a harmonic potential and the side chain atoms were unrestrained. The 17 NMR solution structures (Krezel et al., 1995) were clustered with the program NMRCLUST (Kelley et al., 1996). Four clusters were found, and a representative structure from the cluster that had 12 members (Fig. 1 a) was chosen. An internal harmonic restraint based on the root-mean-square deviation (rmsd) relative to the experimental NMR structure (Krezel et al., 1995) of AgTx2 was applied to the backbone, and the side chains of the toxin were unrestrained. This rmsd energy restraint allows complete freedom of translation and rotation while maintaining the conformation of the backbone of the toxin close to the reference structure. The strength of the rmsd backbone restraints was decreased from 100 to 5 kcal/mol/Ų during the docking. All the docking trajectories were started with different initial conditions. The toxin was given a completely random orientation and was positioned randomly in a region of 5 by 5 by 1 Å centered along the pore axis at a distance of 20 Å from the channel on the extracellular side. Hydrogen atoms were not included, and electrostatic interactions were ignored in the early stage of the docking simulation. Atomic displacements along a fictitious fourth dimension were allowed to reduce the core-core repulsion between the toxin and channel, thus enabling easier structural rearrangements of the side chains at the protein-protein interface (i.e., this MD in the fourth dimension literally allows atoms go through each other in 3D; van Schaik et al., 1993). All atoms were included, and electrostatic interactions were taken into account toward the end of the docking simulations. The initial temperature was 500 K. The temperature and the position of AgTx2 along the fourth dimension were gradually decreased to generate structures of the complex at 300 K in 3D. The time step was 1 fs, and the total length of one docking simulation was 10 ps. The distance restraint force constant, K, was gradually increased from 10, 10, and 10 kcal/mol/Ų to 60, 40, and 15 kcal/mol/Ų for the strong (s), medium (m), and weak (w) restraints, respectively, during the docking simulations (Table 1). Models with distance restraints energies greater than 100 kcal/mol, corresponding to a deviation of ~1 Å on the distance restraints, were rejected. The best models resulting from the docking trajectories (in terms of distance restraints energies) were further refined using a simulated annealing protocol with electrostatics during which the temperature was decreased from 1000 to 300 K, and the fourth coordinate of AgTx2 was decreased to zero. To relax the protein-protein complex, the distance restraints were gradually weakened as the temperature decreased. Hydrogen atoms were restrained with the SHAKE algorithm (Ryckaert et al., 1977), and the time step was 2 fs. The total length of one simulated annealing was 28 ps. All the MD calculations were performed using the CHARMM program version c28a3 (Brooks et al., 1983).

To test of the docking simulation protocol, Barstar was docked onto Barnase according to the same protocol that was used for determining the AgTx2-Shaker complex. The Barnase-Barstar complex is very suitable as a test system because a high-resolution x-ray structure is available (Buckle et al., 1994; Guillet et al., 1993) as well as binding free energies from several single and double mutants. The set of distance restraints used is given in Table 2. Each docking simulation was started with a randomized position of Barstar along the principal axis above the interface and a random orientation of Barstar relative to Barnase. The Barnase backbone atoms were restrained relative to the crystal structure using a harmonic potential. An internal harmonic restraint based on the rmsd relative to the crystal structure of the complex was applied to the Barstar backbone. The side chains of both Barnase and Barstar were fully flexible. Fifty docking simulations were performed using the described protocol and models with a distance restraints energy <50 kcal/mol (35 of 50) were refined with the simulated annealing protocol. The rmsd of the backbone of Barstar relative to the crystal structure of the complex was 0.7 Å for the best model. Most side chain conformations at the Barnase-Barstar interface in the resulting model are also in good agreement with the x-ray structure with an rmsd of 1.8 Å for all heavy atoms. Thus, given that there are enough distance restraints between the two proteins, the docking protocol is expected to give satisfactory results.

Thermodynamic double mutant cycles

In a thermodynamic double mutant cycle, the interaction G_int between two residues "a" and "b," respectively, in protein "A" and "B" is determined by considering the difference between the single and double mutant relative binding free energy (Ranganathan et al., 1996; Schreiber and Fersht, 1995),

(2)

in which G_bind (a, b), G_bind(a*, b), G_bind (a, b*), and G_bind (a*, b*) represent the binding free energy of the wild-type, single, and double mutant proteins.  is a measure of the nonadditivity of the mutations "a  a*" and "b  b*" on the binding affinity (Ranganathan et al., 1996; Schreiber and Fersht, 1995). If the interactions other than those between a residue pair canceled out, the coupling constant would be roughly related to the difference in direct interaction energy (Horovitz et al., 1990; Serraro et al., 1990).

A continuum solvent model provides a useful and computationally inexpensive approach for calculating the binding free energies of the various mutants (Norel et al., 2001; Novotny et al., 1997; Olson, 1998; Olson and Reinke, 2000; Sharp, 1998). To make progress, the binding free energy G_bind is expressed as

(3)

in which G_np and G_elec are the nonpolar and electrostatic contributions, respectively. The electrostatic contribution is calculated from the Poisson-Boltzmann equation. The nonpolar contribution to the binding free energy is empirically written as a fraction of the van der Waals interactions E_vdW upon formation of the complex between protein "A " and "B."

The value of , as well as the dielectric constant assigned to the protein interior, _prot, was optimized empirically using the experimental data about the 13 single mutations of the Barnase-Barstar complex (Ranganathan et al., 1996; Schreiber and Fersht, 1995). For each mutation, the structure was relaxed using a short MD simulation followed by an energy minimization for all atoms at a distance less than 6 Å from the mutated residue(s), while the rest of the complex was kept fixed. The Poisson-Boltzmann equation was solved numerically using the PBEQ module (Im et al., 1998) implemented in the program CHARMM (Brooks et al., 1983). A set of atomic Born radii, calibrated and optimized to reproduce the electrostatic free energy of the 20 amino acids in MD simulations with explicit water molecules, was used (Nina et al., 1997). An excellent agreement was obtained with _prot = 12 and  = 0.17; the rmsd between calculated and measured binding free energies, , was equal to 1.3 kcal/mol. Incorporating the influence of the solvent accessible surface area (SASA) did not improve the results. The small value of  reflects the fact that the short-range dispersive interactions in the complex are partly (but not completely) compensated by protein-solvent interactions in the dissociated state. The relative binding free energy of the single mutants are shown in Fig. 2 (top). The G_bind = [G_bind (mut)  G_bind (wt)] for 14 double mutations in the x-ray structure yields a  value of 1.5 kcal/mol. The G_int for the double mutants are shown in Fig. 2 (bottom). The agreement with experimental values is excellent with a  of 1.4 kcal/mol.

View larger version (45K): [in this window] [in a new window]   FIGURE 2   (Top) Gbind (kcal/mol) for 13 single mutations in Barnase-Barstar. (Bottom) Gint (kcal/mol) for 14 Barnase-Barstar residue pairs. Experimental values (Schreiber and Fersht, 1995) are shown as black bars and the values calculated from the x-ray structure (Buckle et al., 1994; Guillet et al., 1993) as light shaded bars.

As a blind test of the applicability of the present continuum solvation model, the binding free energy of the lysozyme dimeric antibody D1.3 complex was examined (Dall'Acqua et al., 1998). As in the case of Barnase-Barstar, the G_bind, calculated for 13 single mutations in either lysozyme or in the antibody were in good agreement with experimental results with a  of 1.0 kcal/mol. Finally, calculations of G_bind and G_int for eight double mutants in the complex resulted in  values of 1.2 and 0.9 kcal/mol, respectively, suggesting that the empirical parametrization is valid for other protein-protein complexes. Such results are comparable with previous calculations by others using similar continuum solvent approximation (Novotny et al., 1997; Olson, 1998; Olson and Reinke, 2000; Sharp, 1998).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Docking and simulated annealing

The docking simulations to determine the structure of the AgTx2-Shaker complex were staged in four sets of calculations. The number of distance restraints was increased progressively in the successive set of calculations, with one restraint in set I, two in set II, seven in set III, and nine in set IV, each set helping to better resolve the structure of the complex by incorporating more restraints. The main results from each series of calculations are summarized in Table 3 and in Fig. 3.

View larger version (42K): [in this window] [in a new window]   FIGURE 3   Distribution of modeled structures for docking set II, III, and IV. Mode I is red, mode I' is green, mode II is yellow, and mode II' is blue in the right figure. The labeling of the monomers in Shaker is shown in the left figure.

A first set of calculations was performed using a single distance restraint involving the side chain of Lys-27. There are several experimental indications that this residue is particularly important for toxin binding. The binding affinity of CTX (a close homolog of AgTx2) decreases threefold upon charge neutralization of Lys-27 and removes the toxin's characteristic voltage dependence in dissociation rate (Goldstein and Miller, 1993; Park and Miller, 1992). Furthermore, ion entry from the internal solution accelerates CTX dissociation from the extracellular face of the channel, suggesting that it competes with K⁺ for a binding site located deeply into the pore (MacKinnon and Miller, 1988; Park and Miller, 1992). These observations suggest that the positively charged amino group of Lys-27 is repelled by K⁺ ions and is thus probably positioned in the pore of the channel (Goldstein and Miller, 1993; Park and Miller, 1992). Lastly, double mutant cycle analysis shows that the coupling constant for the Lys-27Met/Tyr-445-Phe is relatively large (Table 1), further indicating that Lys-27 could bind directly into the pore (Ranganathan et al., 1996).

Although plausible, it remains nevertheless uncertain whether the molecular shape of the two proteins can sterically accommodate the binding of Lys-27 in close proximity to the pore. Can the lysine side chain actually reach into the pore or not? To examine this possibility, 200 models of the toxin-channel complex were generated following the simulation protocol described in Theory and Methods. A single distance restraint between the positively charged amino group of the Lys-27 side chain and the carbonyl oxygens of Tyr-445 was applied, positioning the positively charged end of the lysine side chain near the outer binding site identified in the x-ray structure (Doyle et al., 1998). This position is also referred to as the S₁ binding site (Bernèche and Roux, 2001). In all of the resulting models, it was observed that the side chain of Lys-27 is located in the pore near the outer binding site, demonstrating that such configurations are sterically possible. Nonetheless, the orientation of the toxin is widely distributed around the pore axis with a slightly higher occurrence observed for the complexes in which Arg-24 is near the residue Asp-431 from one of the four channel subunits. Although these initial calculations demonstrate that the docking procedure is able to explore extensively over all orientations of the toxin bound at the extracellular surface of the channel, they show that the restraint on Lys-27 alone is not sufficient for determining an accurate model. Presumably, the complementarity between the molecular surface of the channel and the toxin is not sufficient to impose a unique orientation for the protein-protein complex. To increase the resolution of the resulting models, it is necessary to introduce additional distance restraints.

A second set of calculations (set II) was performed, including an additional restraint between Arg-24 and Asp-431. This pair of residue exhibits the largest coupling observed experimentally in the double mutant experiments (ln() = 7.3, see Table 1), probably reflecting the formation of a salt bridge in the toxin-channel complex (Hidalgo and MacKinnon, 1995). Consequently, a restraint imposing a maximum distance of 3.0 Å between any pair of oppositely charged atoms was applied (i.e., N₁ or N₂ of Arg-24 with O₁ or O₂ of Asp-431). In principle, Arg-24 can bind alternatively with Asp-431 on any of the four monomers of Shaker, although the symmetrically related complexes are structurally equivalent and provide no additional information. To avoid the unnecessary complication of the fourfold symmetry, the distance restraint between Arg-24 and Asp-431 was assigned only to monomer A of Shaker; 100 models were generated following the simulation protocol. A total of 80 models satisfying the distance restraints were kept, and the remaining models were rejected. For this set of 80 structures, the rmsd of the toxin backbone relative to the average structure is 3.7 Å. The resulting models allowed the assignment of a few additional unambiguous distance restraints corresponding to the strongest coupling in the double mutant data. In particular, it was observed that Phe-25 is within 4 Å of Met-448 of monomer A and that Pro-37 is less than 4 Å from Val-451 of monomer B.

Even though the experimentally observed coupling between Ser-11 and Asp-431 is relatively large, these residues were in close proximity in none of the generated models. In fact, further analysis indicated that satisfying both the Arg-24-Asp-431 and Ser-11-Asp-431 distance restraints was not possible even though they both exhibit strong coupling in the double mutant data; 50 models were generated with the Lys-27-Tyr-445 distance restraint and one additional ambiguous restraint between Asp-431 and either Ser-11 or Arg-24 (the residue that happened to be nearest). On the basis of these models, it was found that Ser-11 was as far as 12 Å away from the nearest Asp-431 when Arg-24 was near Asp-431. Conversely, when Ser-11 was close to Asp-431, Arg-24 was 8 to 15 Å away from the nearest Asp-431. This computational experiment suggests that the observed Ser-11/Asp-431 coupling does not reflect a short distance in this case. For this reason, no distance restraint was applied between Ser-11 and Asp-431 in the docking simulations.

On the basis of the analysis of the previous models, 150 additional models were generated with four monomer-specific restraints and three additional ambiguous restraints (set III). A total of 109 models satisfying the distances restraints were kept, and the remaining models were rejected from the analysis. The rmsd of the toxin backbone for the ensemble of models relative to the average structure is 3.0 Å, which is slightly smaller than that of the previous set of calculations. The resulting models allowed the assignments of additional monomer-specific contacts between the toxin the channel: Gly-10 is close to Thr-449 of monomer D, and close to Phe-425 in either monomer C or D, whereas Ser-11 is close to Asp-447 in either monomer C or D. The models were also analyzed to identify additional contacts detected in the double mutants data experiments, which could contribute to better define the structure of the toxin-channel complex. It was found that Phe-25 is close to Val-451 in monomer A, whereas Met-29 is close to Thr-449 in either monomer C or D. Couplings on the order of 0.5 k_BT were measured in the double mutant cycles for these pairs of residues (Ranganathan et al., 1996).

To better resolve the resulting protein-protein complexes, two separate sets of 50 models were generated using seven monomer-specific restraints and two ambiguous restraints (set IV): a first set in which Ser-11 was restrained to Asp-447 in monomer C and a second set in which it was restrained to Asp-447 in monomer D. The rmsd of the toxin backbone for the resulting 100 models is only 1.8 Å. From these models it was possible to distinguish four different binding modes: I, I', II, and II'. The rmsd for the individual mode is only 0.5 Å, which indicate that they are distributed quite narrowly. The four binding modes occurred with similar frequencies and with similar distance restraint energies and, thus, appear to be equally plausible. The four binding modes can clearly be observed in Fig. 3 (right). In particular, mode I (red) and mode I' (green) are quite similar to each other, and similarly for mode II (yellow) and mode II' (blue). However, it can be observed that the orientation of the toxin differ by a rotation on the order of ~20° around the z axis of mode I (or I') relative to mode II (or II'). The relative rmsd is on the order of 1.2 Å. Further structural refinement using simulated annealing with all interactions on the eight best structures of each binding mode (lowest restraint energies) decreased the backbone rmsd of AgTx2 relative to the average structure from each of the four binding modes to 0.4 Å. The significant structural difference between the four binding modes indicates that the procedure has reached its intrinsic limit. Further discrimination between the different models is thus needed to converge toward a valid model of the AgTx2-Shaker complex.

Unrestrained MD simulations with explicit solvent

The stability of the four binding modes was examined by generating 3 ns MD trajectories of the toxin-channel complex solvated by explicit water molecules. The strengths of the distance restraints were gradually reduced to zero at 0.4 ns; internal backbone restraint on AgTx2 were gradually reduced from 1.6 to 1.8 ns. The time evolution of the AgTx2 rmsd from the starting structure and the minimum distances between the restrained residue pairs are shown in Fig. 4. After 2.2 ns the protein complexes of mode I and II were both stable and nondrifting with a backbone rsmd from the initial structure around 1.6 Å. All distances between the initially restrained residue pairs were also stable and nondrifting (Fig. 4, left). In contrast, mode I' and II' appeared to be much less stable. In mode I' (Fig. 4, top, right), the AgTx2 backbone rmsd is slightly higher (2.0 Å) and the N-terminal part of the helix (around residue 10 to 12) is opening up toward the end of the simulation. In the case of mode II', there is an abrupt increase in the AgTx2 backbone rmsd when the internal restraint is removed (Fig. 4, bottom, right) and the rmsd is drifting throughout the simulation. The G10-F425 and Ser-11-D447 distances in mode I' and II' are significantly larger than the corresponding distances in mode I and II.

View larger version (46K): [in this window] [in a new window]   FIGURE 4   MD simulations of one annealed model from each binding mode of the toxin-channel complex with explicit solvent molecules. Each colored curve shows the time evolution of the minimum distance between the initially restrained pairs in AgtTx2-Shaker for the four binding modes. The violet curve is the time evolution of the AgTx2 backbone rmsd from the initial structure, and the other curves are minimum distances between residue pairs (as defined in Table 1). The simulated systems were solvated with a 25-Å radius sphere of water molecules centered at Tyr-445. The strength of the restraints was gradually decreased until the molecular complex was freely simulated at 0.4 ns. The nonbonded interactions were truncated at 10 Å; the length of bonds involving hydrogen atoms were kept fixed using the SHAKE algorithm (Ryckaert et al., 1977), and the time step was 2 fs.

At the end of the simulations the AgTx2 backbone rmsd relative to the NMR structure are on the order of 1 to 1.5 Å. The small differences are mostly located at the turn between residues 29 to 31. Examination of the toxin-channel interface shows that all side chains have unique conformations within one binding mode, except those at the protein-solvent interface. The side chain conformations of AgTx2 are generally the same as in the NMR structure with Lys-27 and Phe-25 being two notable exceptions. In all the AgTx2-Shaker complexes, the side chain of Lys-27 adopts an extended conformation, pointing directly into the pore. However, such a conformation is observed in none of the NMR structures (Krezel et al., 1995). Similarly, the rotameric state of the Phe-25 side chain differs from the one observed in the NMR structure; the ₁ and ₂ dihedral angles of the Phe-25 side chain changed from 178° and 66°, respectively (in the NMR structure), to 159° and 80° in mode I and to 84° and 60°, respectively, in mode II. These conformational changes of AgTx2 upon docking to Shaker show that there is partly an induced fit in toxin-channel interactions. These observations highlight the importance of allowing side chain flexibility in searching for the optimal conformation of protein-protein complexes, an aspect that is generally neglected in docking algorithms for the sake of computational efficiency (for example, see Camacho and Vajda, 2001; Palma et al., 2000; Ritchie and Kemp, 2000). Side chain flexibility of toxins has been ignored in previous rigid-body BD docking studies (Cui et al., 2001; Cui et al., 2002).

The magnitude of the SASA buried upon formation of the AgTx2-Shaker complex is roughly 1000 Å², a value that is comparable with that of other known protein-protein complexes (Glaser et al., 2001; Jones and Thornton, 1996). The SASA covered by the toxin in the case of mode I and I' is roughly 350 Å² for monomer A, 250 Å² for monomer B, 250 Å² for monomer C, and 150 Å² for monomer D. In the case of mode II and II', the SASA is 360 Å² for monomer A, 230 Å² for monomer B, 120 Å² for monomer C, and 240 Å² for monomer D. The variations in the SASA reflect the rotation of the toxin in the different binding modes.

Binding free energies and coupling constants

To further assess the validity of the toxin-channel complexes, the binding free energies of the single and double mutants were calculated using a continuum solvation approximation. The total electrostatic and nonpolar contributions to the binding free energy (G_elec + G_np) of the four binding modes were calculated as an average over eight configurations, taken near the end of the MD simulations. In all cases, G_elec + G_np was very similar for each binding mode, 38 ± 2 kcal/mol. One may note that such value is roughly consistent with the experimental dissociation constant in the nanomolar range (Ranganathan et al., 1996) once the loss of entropy caused by the translational and orientational restriction (approximately +7 kcal/mol) as well as the immobilization of side chains with rotatable bonds at the toxin-channel interface (approximately +16 kcal/mol; Doig and Sternberg, 1995) have been accounted for. However, the calculated absolute free energies are much too similar to allow discrimination between the different binding modes.

To further discriminate among the binding modes, we calculated the difference in binding free energies between mutated and wild-type complexes (G_bind) for 18 single mutations in each binding mode. The G_bind were calculated from an average over four configurations, taken near the end of the MD simulations; the standard deviations are on the order of ±0.2 kcal/mol. The mutations involving the four tyrosines of the selectivity filter (Y445F) were not included in those calculations because their influence on toxin binding is expected to be mediated by subtle and complex effects involving the carbonyl oxygens that would be inaccurately modeled by the continuum solvent model. The calculated and measured G_bind for the four binding modes are shown in Fig. 5. Overall, the agreement with experiments is very reasonable for all modes. The results are slightly better results in the case of mutations in AgTx2 than in Shaker, reflecting the fact that a mutation in Shaker involves one point mutation in each one of the four subunits, which can lead to an accumulation of small errors. Because the coupling constants  represent the raw information extracted from experimental double mutants cycles used in the docking simulations, they are of particular interest in trying to assess the validity of the models. Errors in the calculated binding free energies of single and double mutants may cancel in the calculations of the coupling constants. Seventeen coupling constants were calculated for double mutants according to Eq. 2. The calculated coupling constants for the four binding modes are shown in Fig. 6. The agreement of the coupling constants with experimental values is significantly better for mode I and II than for mode I' and II' with rmsd values of 1.6 (0.9 kcal/mol) and 1.3 (0.8 kcal/mol), respectively.

View larger version (45K): [in this window] [in a new window]   FIGURE 5   Calculated and experimental Gbind for 18 single mutations in either AgTx2 or Shaker.

View larger version (50K): [in this window] [in a new window]   FIGURE 6   Calculated and experimental coupling constants, ln(), for 17 double mutants. No distance restraints were imposed on the residue pairs that are labeled with an asterisk.

Generally the magnitude of ln() is much smaller for mode I' and II' than for mode I and II. The low values of the coupling constants are especially pronounced for mode II' (blue bars in Fig. 6). Two clear examples are the ln() values for Asp-447-Glu/Ser-11-Ala and Phe-425-Gly/Gly-10-Val, which are close to zero in both modes I' and II', whereas they have values close to the experimental in the two other modes. This reflects the fact that the toxin is slowly shifting away from the surface of Shaker, consistent with the observation that the MD simulations of both mode I' and II' are less stable than the two other modes. All these results suggest that mode I' and II' can be ruled out as candidates for the correct binding mode.

It is particularly encouraging that the values of ln() for residue pairs where distances between the residues were not restrained during the MD simulations (marked with an asterisk in Fig. 6) are also in good agreement with the experimental values for modes I and II. The Asp-431/S11 pair is particularly interesting. Even though a distance restraint could not be satisfied in the docking simulations (see above), the calculated value of ln() have the correct sign (but are underestimated). Overall, the agreement with the experimental data is slightly better for mode II than mode I. This is mostly due to the value of ln() for Met-448-Ala/Phe-25-Ala. In mode II, these mutations give rise to a ln() close to the experimental, whereas in mode I it is close to zero. Although the Met-448-Phe-25 distance is roughly 4.5 Å for both binding modes (blue curves in Fig. 4), the phenyl ring of Phe-25 is oriented differently in the two binding modes (see above). The plane of the phenyl ring is perpendicular to the pore axis is mode I, whereas it is more parallel in mode II. The better agreement for mode II suggests that the orientation of the Phe-25 phenyl ring should be almost parallel to the pore axis rather than perpendicular.

All of the AgTx2-Shaker residue pairs that are less than 4 Å apart from each other in either one of the two binding modes are listed in Table 4. The particular interresidues distances, which differ markedly between mode I and mode II, are highlighted in bold. The two regions of AgTx2 that differ most between the binding modes are the N-terminal part of the helix (residues 9-12) and the -turn between residues 28 and 32. For example, the Thr-9-Phe-425 and Ser-11-Thr-449 distances are both ~3 Å larger in mode II than in mode I (Table 4). The weak couplings that were measured between Thr-9-Phe-425 and Ser-11-Thr-449 should reflect relatively large interresidue distances. In both cases, the distance between the pair of residue is smaller in mode I, suggesting that mode II is more consistent with the experimental data. Conversely, the distance between Asn-30-Phe-425 is almost 5 Å larger in mode I than in mode II. The relatively strong coupling between this residue pair (ln() = 1.1) is in much better accordance with the distance between the residues as in mode II. Fig. 7 shows the AgTx2-Shaker complex in binding mode II and summarizes some important contacts found between the toxin and the channel.

View larger version (61K): [in this window] [in a new window]   FIGURE 7   (Top) Summary of important contacts between AgTx2 and Shaker in binding mode II. (Bottom) Two views of the AgTx2-Shaker complex (mode II) rotated 90° around the pore axis relative to each other. The rear and front monomers have been removed for easier visualization and the letters refer to Shaker monomers. The potassium ion in the selectivity filter is shown as a magenta sphere, and the side chains have the same colors as in the figure above; (left) Lys-27-Tyr-445-O, red; Arg-24-Asp-431(A), green; Phe-25-Met-448(A), magenta; Lys-19-Glu-422(A), orange; and Asn-30-Phe-425(C), cyan. (Right) Lys-27-Tyr-445-O, red; Pro-37-Val-451(B) and Thr-9-Glu-422(D) blue; Lys-38-Glu-422(B), magenta; Ser-11-Asp-447(D) yellow; and Pro-12-Phe-425(D), cyan.

For some specific residue pairs there are obvious problems with the calculated couplings. In particular, the calculated coupling constants involving the Phe-425-Gly or the Thr-449-Cys mutations have the wrong sign for all binding modes. Although the reasons for such discrepancies are unclear, it might reflect the intrinsic limitations of the simple continuum solvent treatment used to calculate the binding free energies, or also indicate that some of the side chains of Shaker are in a wrong conformation (e.g., in particular the bulky aromatic ring of Phe-425). Difficulties arising from the 3D model of Shaker should be expected. Despite the fact that the sequence similarity of Shaker and KcsA is generally very high (with no insertion or deletion), the residues in the turret region (gray in Fig. 1 b) have a very low sequence similarity (two of eight residues). For this reason, it is likely that these regions of the 3D model are not represented accurately. For example, in mode II, the charged groups of Lys-16 and Lys-19 form strong contacts with Glu-422 and Ser-424, respectively. In mode I, the conformations of these side chains are somewhat different as Lys-16 points out in the solution and Lys-19 instead interacts with the backbone carbonyl of Gly-422. However, the experimental values of ln()for Lys-16-Glu-422 and Lys-38-Glu-422 are only 0.5 and 0.2, respectively (Hidalgo and MacKinnon, 1995), which suggests that no salt bridge is formed between these residues. Lastly, the G_bind for the Asn-30-Ala mutation is underestimated by 3 to 4 kcal/mol in all four binding modes. No distance restraints were available to determine the orientation of this polar side chain, although its large effect on the association constant indicates that it should interact strongly with the extracellular surface of Shaker. In the toxin-channel complex Asn-30 is relatively close to Asp-431. To examine the occurrence of an interaction, a short simulated annealing was performed with an additional distance restraint between Asp-431-O_1/2 and Asn-30-N₂. Although such an additional restraint could not be satisfied without significant perturbation of the complex in the case of mode I, it could easily be satisfied in the case of mode II. In this case, the backbone of AgTx2 backbone moved only slightly (rmsd relative to the previous structure is 1.6 Å), the changes affecting mostly the loop region between residues 28 to 32. In the resulting structure, the Asn-30 side chain is 3.2 Å away from Asp-431 and also within 5 Å from Asp-447. This computational experiment shows that Asn-30 could interact more strongly with Shaker without any significant structural changes of the complex in the case of mode II but not in the case of mode I.

Importantly, it is possible to use the present atomic models to guide the design of experiments aimed at further refining the structure of the AgTx2-Shaker complex. In particular, Pro-12 is located at the N terminus of the helix in AgTx2 where the structural difference between the binding modes is the largest. The residue forms various contacts with monomer D of Shaker (Table 4), and the difference in distance between the binding modes for a given residue pair is ~3 Å. The different contacts in modes I and II should give rise to different couplings and help to discriminate between the modes. Hypothetically, stronger couplings should be found for Pro-12-Ser-424, Pro-12-Met-448, and Pro-12-Thr-449 in the case of mode I, whereas a stronger coupling should instead be found for the pair Pro-12-Phe-425 in the case of mode II. Similarly, Gly-10 is in contact with Asp-447 in mode I, whereas is it 6.5 Å away from the closest Asp-447 in mode II. Thus, a strong coupling would be expected for Gly-10-Asp-447 only in the case of mode I. Finally, it would be interesting to determine the coupling of Asn-30 with Asp-431 and with Glu-422.

Although the validity of mode II appears to be better than mode I, it is somewhat surprising that the calculated properties are in such good agreement for two significantly different AgTx2-Shaker complexes. Both mode I and II are stable during the 3-ns MD simulations with explicit water molecules and have very similar calculated coupling constants (slightly better for mode II). In view of these results, it is tempting to speculate that AgTx2 might not bind to Shaker in a unique way. From a functional point of view, the existence of such multiple binding modes might effectively increase the capture radius of AgTx2 for the surface of Shaker, which would result in a more efficient toxin. Although we see no compelling reasons to reject the possibility of multiple binding modes, it seems more likely that the lack of convergence toward a unique structure reflects the intrinsic limitations of the set of coupling constants that is currently available. Ultimately, additional experiments are required to better resolve the structure of the AgTx2-Shaker complex. The usefulness of the present models will be in guiding the design of future experiments.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The purpose of this work was to determine the structure of the neurotoxin AgTx2 in complex with the pore of Shaker potassium channel. A docking protocol using distance restraints between AgTx2 and Shaker derived from experimental double mutant data was applied for the structure determination. A series of dockings simulations was performed during which the number of distance restraints was successive increased. The series converged into four different complexes of AgTx2-Shaker, corresponding to two similar binding modes (I, I' and II, II', respectively). The orientation of the toxin in the two binding modes differ quite significantly by a rotation of ~20° around an axis parallel to the pore.

To discriminate between the binding modes, 3-ns MD simulations of the complex were performed with explicit water molecules. The simulations of mode I and II were both stable, with a nondrifting AgTx2, whereas mode I' and II' were significantly less stable. Thereafter, binding free energies of 18 single mutant complexes relative to the wild-type complex (G_bind) were calculated from structures taken near the end of the MD simulations. The agreement with experimental data was reasonable, but no significant difference was observed between the four modes. The subsequent calculation of coupling constants () from 17 double mutants turned out to be more helpful for discriminating between the binding modes. The agreement with experimental data was much better in the case of mode I and II relative to I' and II', suggesting these two modes do not represent good models of toxin-channel complex. In both models, the side chain of Lys-27 binds at the outer binding site for K⁺ and plugs the pore. Further analysis suggests that toxin-channel contacts for some residues located in the regions where the structural difference between the two modes is largest (e.g., Thr-9, Ser11, and Asn-30) are more consistent with the experimentally measured coupling constants in the case of mode II, although the lack of clear convergence toward a unique structure indicates that additional coupling constants are required to better resolve the structure of the AgTx2-Shaker complex. The present modeling suggests that those involving Pro-12 and Asn-30 on the toxin and Ser-424, Met-448, Thr-449, and Phe-425 on Shaker would be particularly informative.

Lastly, the present study illustrates well the power of double mutant cycles in combination with computational methods for the structural characterization of macromolecular complexes. This strategy may be helpful when a high-resolution structure of a complex cannot be obtained by NMR or x-ray crystallography even though the structure of the constituents is known.