INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION CONCLUSION REFERENCES

The determination of the three-dimensional structure of the KcsA K⁺ channel of Streptomyces lividans at an atomic resolution of 3.2 Å, using x-ray crystallography (Doyle et al., 1998), represents an extraordinary opportunity for understanding ion permeation in biological membrane channels at the molecular level. The KcsA K⁺ channel shows sequence similarity to all known K⁺ channels (Doyle et al., 1998; Schrempf et al., 1995; Cortes and Perozo, 1997). In particular, the amino acid sequence of the central domain is very similar to the segment S5, H5, and S6, which is conserved in eukaryotic voltage-activated channels such as the Shaker of Drosophila melanogaster (Doyle et al., 1998). In the pore region, the channel exhibit the typical "signature sequence," TVGYG, common to all K⁺ channels (Heginbotham et al., 1992, 1994). Its permeation properties are remarkably similar to those of standard K⁺ channels (Schrempf et al., 1995; Cuello et al., 1998; Heginbotham et al., 1998, 1999). Furthermore, based on a combination of structural and functional data on scorpion neurotoxin, it can be concluded that the KcsA K⁺ channel is structurally very similar to eukaryotic K⁺ channels (MacKinnon et al., 1998). Thus investigations of KcsA are expected to help in the understanding of a large class of biologically important channels.

The general architecture of KcsA, shown in Fig. 1, reveals the basic mode of operation of K⁺ channels. The structure is made of four identical subunits of three -helices dispersed symmetrically around a common axis corresponding to the pore. The narrowest part of the pore, formed by the backbone carbonyl oxygens of Thr⁷⁵-Val⁷⁶-Gly⁷⁷-Tyr⁷⁸, acts as a selectivity filter for K⁺ ions. It is only 12 Å long and is located near the extracellular side. Two ion binding sites separated by ~7.5 Å are observed in the selectivity filter, the outer site, on the extracellular side and the inner site, on the intracellular side. At the level of the bilayer center, the radius of the pore increases to form a 5-Å-radius water-filled cavity lined by nonpolar residues. The widening of the pore contributes to the maintenance of a high throughput rate of ions by minimizing the distance over which the ions interact strongly with the channel. Four -helices of 13 residues, the "pore helices," are positioned at an angle of ~45° with respect to the pore axis, so that their COO terminus points toward the center of the cavity. Detailed calculations based on the Poisson-Boltzmann equation show that the pore helices electrostatically stabilize a monovalent cation in the cavity center (Roux and MacKinnon, 1999). The cavity and pore helices thus help to overcome the high electrostatic energy barrier presented by the low dielectric membrane (Parsegian, 1969) to achieve a large throughput rate.

View larger version (46K): [in this window] [in a new window]   FIGURE 1   Schematic view of the KcsA K+ channel. The extracellular side is at the top and the intracellular side is at the bottom. The main structural elements are the outer helix (residues 27-51, in red), the pore helix (residues 62-74, in yellow), the selectivity filter (residues 75-79, in green), and the inner helix (residues 86-112, in blue). Three K+ ions, located at the outer and inner sites of the selectivity filter and at the center of the central cavity, are also shown (purple).

Although the crystallographic structure does provide essential information about the KcsA channel as well as the ions and water in the pore, further considerations are necessary to understand the mechanism governing ion permeation at the microscopic level. Because proteins are very complex dynamic macromolecular systems, a static picture of the three-dimensional spatial arrangement of the atoms is not sufficient to understand how they function (Brooks et al., 1988). The difficulties are particularly acute in the case of a membrane transport protein such as an ion channel, because the underlying microscopic processes involve energetic and dynamic factors that are not directly accessible to experiment (Hille, 1992). In particular, little or no information is available about the transient configurations of the channel, ions, and water molecules occurring during the permeation process. Furthermore, some details of the present three-dimensional structure remain uncertain because of the moderate resolution of the x-ray data. In particular, the side chain of Glu⁷¹, an important residue in the vicinity of the selectivity filter, could not be completely resolved in the electron density (Doyle et al., 1998). Its atomic coordinates are undetermined and its protonation state is not known. Last, alternative configurations of the K⁺ ions and water molecules in the selectivity filter, deduced experimentally by analyzing difference Fourier maps of crystals containing Rb⁺ and Cs⁺ relative to those with K⁺, are possible on the basis of the current data (Doyle et al., 1998).

The goal of this paper is to begin to explore the structure and dynamics of the KcsA K⁺ channel with molecular dynamics (MD) simulations based on detailed atomic models. Our aim with these calculations is to complement the information that is currently available from experiments and, ultimately, to refine in our understanding of K⁺ channels. In particular, the simulations are used to examine the configuration of the side chain of Glu⁷¹ and of the K⁺ ions and water molecules in the pore, for which the information is limited because of the moderate resolution of the x-ray data. In addition, the microscopic fluctuations of the channel as well as other factors that cannot easily be accessed experimentally (but are important for an understanding of the function of KcsA) are characterized.

Current MD methodology has reached the point where one can generate trajectories of realistic atomic models of complex biological systems. In recent years, the approach has provided a great wealth of information about transmembrane ion channels, such as gramicidin (Woolf and Roux, 1996; Chiu et al., 1999), LS3 (Zhong et al., 1998), OmpF Escherichia coli porin (Tieleman and Berendsen, 1998), alamethicin (Tieleman et al., 1999), and KcsA (Guidoni et al., 1999; Shrivastava and Sansom, 1999). In the present paper, we have generated dynamical trajectories of detailed atomic models of the KcsA K⁺ channel embedded in a fully solvated bilayer membrane with explicit lipids and aqueous salt solution.

In the next section, the atomic model and the computational methodology are described in detail. The results are then analyzed and discussed. The paper concludes with a brief summary of the main results.

	THEORY AND METHODS

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION CONCLUSION REFERENCES

Microscopic system, potential function, and simulation procedure

The simulation system represents an atomic model of the KcsA channel embedded in dipalmitoylphosphatidylcholine (DPPC) surrounded by a 150 mM KCl aqueous salt solution. The model comprises KcsA (four submits of 97 amino acids for a total of 5292 atoms), 112 DPPC (46 and 66 in the top and bottom layers, respectively), 6532 water molecules, 3 K⁺ in the pore (two in the selectivity filter and one at the cavity center), and 12 K⁺ and 23 Cl in the bulk solution. The total number of atoms in the system is slightly above 40,000. The entire system is electrically neutral. The membrane normal is oriented along the Z axis, and the center of the bilayer is at Z = 0. KcsA was oriented with the pore along the Z axis and the selectivity filter located near the upper layer of the membrane.

The calculations were performed using the academic version c27a1 of the biomolecular simulation program CHARMM (Brooks et al., 1983). The all-atom potential energy function PARAM-22 for protein (MacKerell et al., 1998) and phospholipids (Schlenkrich et al., 1996) was used. The TIP3P potential was used for the water molecules (Jorgensen et al., 1983). The Lennard-Jones (LJ) parameters for K⁺ and Cl were adjusted to yield the experimental solvation free energy in bulk water (Roux, 1996). Periodic rectangular boundary conditions were applied in all directions to simulate a multilayer system of membranes extending in the XY plane. The Z dimension of the unitary cell was allowed to vary according to the constant pressure and temperature thermodynamic ensemble with fixed surface area (CPTA) (Feller et al., 1995, 1997). The dimensions of the periodic system are 72 × 72 × 78 Å³. The electrostatic interactions were computed with no truncation, using the particle mesh Ewald (PME) algorithm with a B-spline order of 4 and a FFT grid of one point per Å (72 × 72 × 80) (Essmann et al., 1995); a real-space Gaussian-width kappa of 0.3 Å¹ was used. The list of nonbonded interactions was truncated at 11 Å, using an atom-based cutoff. The nonbonded van der Waals and real-space electrostatic interactions were smoothly switched off at 8-10 Å. Such a treatment of electrostatic interactions was essential for investigating a molecular system conducting ions. The trajectory was generated in the constant pressure and temperature ensemble (CPT) with extra degrees of freedom for a Langevin thermostat and piston (Feller et al., 1995). The SHAKE algorithm (Rychaert et al., 1977) was used to fix the length of all bonds involving hydrogen atoms, and the equations of motion were integrated with a time step of 2 fs. The coordinates were saved every 0.1 ps.

Construction of the atomic system

KcsA structure

View larger version (40K): [in this window] [in a new window]   FIGURE 2   Detail of the selectivity filter from a single monomer of KcsA. The side chain of Glu71 was modeled to form a strong hydrogen bond with the carboxylate group of Asp80 near the extracellular surface after discussions with R. MacKinnon and co-workers (private communication). The hydrogen bonds that were restrained during the early stage of the equilibration are marked by thick lines (no restraints were applied to the salt bridge between Asp80 and Arg89 during the equilibration).

Ions and water molecules in the pore region

View larger version (12K): [in this window] [in a new window]   FIGURE 3   The four initial configurations of ions and water molecules in the selectivity filter (C1, C2, C3, and C4) that were considered. The inner ion is in the upper site in configurations C1 and C2 and in the lower site in C3 and C4.

DPPC membrane

Aqueous salt solution

Equilibration protocol

The construction and equilibration of the simulation systems are described in Fig. 4. A first system was initially assembled in configuration C1 (see Table 1). The total length of the equilibration period is 600 ps. The system was first thermalized with Langevin dynamics at constant volume. A friction of 3.0 ps¹ was used for the nonhydrogen atoms. The last 200 ps of equilibration was performed with a constant area-isothermal isobaric algorithm with a pressure of 1.0 atm and a temperature of 315 K (Feller et al., 1997).

View larger version (14K): [in this window] [in a new window]   FIGURE 4   Equilibration and simulation protocol for the current calculations. The first 400 ps of equilibration was at constant volume dynamics (NVE), while the remainder of the equilibration and the production was at constant pressure and temperature (NPT). In the first part of the equilibration, harmonic restraints (applied to the center of mass of the polar headgroups, the protein backbone, and the ions) were gradually decreased to allow a smooth relaxation of the system. In the second part of the equilibration the different initial configurations of ions and water molecules in the selectivity filter were constructed and refined. Hydrogen bond restraints were applied to the selectivity filter (see Fig. 2) and gradually decreased.

The purpose of the initial stage of equilibration is to relax the membrane lipids and the aqueous salt solution. For the first 200 ps of equilibration the temperature was 330 K, which is above the gel-liquid phase transition temperature of DPPC (Gennis, 1989). In the remainder of the equilibration, the average temperature of the system was set at 315 K, which is near the gel-liquid phase transition of pure DPPC membrane. Progressively decreasing structural energy restraints were applied to ensure a smooth equilibration of the system (see Fig. 4 for details). As in previous work (Woolf and Roux, 1994, 1996; Bernèche et al., 1998), the displacement of the center of mass of the lipid headgroups away from Z = ± 19 Å was restricted with planar harmonic restraints, and penetration of water molecules in the hydrocarbon region of the bilayer (|Z| < 11 Å) was prevented by using a half-harmonic repulsive planar potential. In addition, harmonic restraining potentials were used to stabilize the hydrogen bonds involving the residues forming the selectivity filter and to limit the displacements of the channel and ions in the pore (see Fig. 2). All of the restraints were gradually reduced to zero at the end of the equilibration period. During the production dynamics, a planar harmonic restraint of 5.0 kcal/mol/Ų was applied to the center of mass of DPPC molecules in the Z direction to prevent any global drift.

A first trajectory of 1 ns at 315 K was generated from the equilibrated C1 (C1/MD-315K). Starting from the same configuration of equilibrated C1, the system was heated to 330 K over a period of 10 ps, and a second trajectory (C1/MD-330K) of 1.5 ns was generated. To examine alternative configurations of the ions and water molecules in the selectivity filter, three additional systems in configurations C2, C3, and C4 (see Table 1 and Fig. 4) were also constructed. Two trajectories of 0.5 ns were generated from C3 at 315 K and 330 K, respectively, using the same procedure (C3/MD-315K and C3/MD-330K). The configurations C2 and C4 were found to be unstable and were not simulated.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION CONCLUSION REFERENCES

The current atomic systems were extensively equilibrated for 600 ps in the presence of progressively decreasing structural energy restraints applied to the ions and the channel. At the end of the equilibration, the systems were completely free of restraints and production trajectories were generated. This extensive equilibration procedure was found to be necessary to obtain meaningful results. Tests simulations with a shorter equilibration exhibited significant distortions relative to the x-ray structure and large displacements of the ions in the selectivity filter. No such large deviations were observed with the extensive equilibration procedure. It is likely that the significant electrostatic forces arising from the 16 carbonyl oxygens pointing toward the two K⁺ ions in the pore are responsible for this behavior.

A total of four trajectories were generated using the initial configurations C1 and C3. During the 1-ns trajectory C1/MD-315K, the system remained roughly near its initial configuration C1, with two K⁺ ions and one water molecule occupying the narrow selectivity filter and one K⁺ ion located in the water-filled cavity. In contrast, a concerted dynamic transition occurred in the pore during the 1.5-ns trajectory C1/MD-330K at 250 ps, after which three K⁺ ions and two water molecules occupied the selectivity filter. A snapshot of C1/MD-330K is shown in Fig. 5.

View larger version (154K): [in this window] [in a new window]   FIGURE 5   Molecular graphics view of the KcsA K+ channel embedded in a DPPC membrane solvated by a 150 mM KCl aqueous salt solution (the K+ ions are purple and the Cl ions are green). This snapshot was taken near the beginning of the trajectory C1/MD-330K.

In the next section, we first analyze the average structural properties of the system based on the trajectories C1/MD-315K and C1/MD-330K. We then analyze the concerted transition in the pore observed in trajectory C1/MD-330K. Finally, alternative configurations of ions and water molecules in the pore are discussed on the basis of trajectories C3/MD-315 and C3/MD-330K.

Average structural properties

During the dynamic trajectories, the structure of KcsA and the DPPC membrane remained very stable. The average root mean square (RMS) deviations of the channel from the crystallographic structure given in Table 2 are relatively modest: the overall RMS deviation of the backbone is 1.9 Å for the whole tetrameric channel, and less than 0.9 Å for the backbone of the main secondary structural elements. The conical shape of the channel is stable. As shown in Table 3, the -helices (outer, pore, and inner) conserve their original orientation with respect to the channel axis. The tilt angle of the pore helix appears to be the most stable, with fluctuations on the order of only 3-5°.

On the intracellular side the inner helices are joined together, and the pore, lined by the hydrophobic residues Thr¹⁰⁷, Ala¹¹¹, and Val¹¹⁵, narrows down to less than 5 Å in diameter. The packing of the inner helices is very stable, and the fluctuations in the diameter of the hydrophobic pore are on the order of 0.6 Å. The central cavity stays filled with ~38 water molecules during the course of the trajectories. No stable contiguous hydrogen-bonded chain of water molecules going from the cavity to the bulk region through the hydrophobic pore is observed. The passage of a cation in this narrow hydrophobic region would require a nearly complete dehydration, which should be highly unfavorable. Nonetheless, such a process has been observed previously in the simulation of Shrivastava and Sansom (1999). Differences in potential function or in simulation methodologies with the current work may be responsible for the different behavior. Shrivastava and Sansom (1999) used the extended atom force field GROMOS (van Giensteron et al., 1999), in which hydrogen atoms are represented as part of the heavy atoms to which they are attached, whereas we used the all-atom PARAM-22 for protein (MacKerell et al., 1998) and phospholipids (Schlenkrich et al., 1996). They truncated the nonbonded electrostatics, whereas we used particle mesh Ewald with no cutoff (Essmann et al., 1995). Truncation of electrostatic interactions, even at a very large cutoff distance, leaves out important ion-ion interactions that could result in simulation artifacts. For example, the distance between the ion at the center of the central cavity and the ion in the outer site is on the order of 17 Å.

Although there are residues with both positively and negatively charged side chains near the mouths of the channel, their spatial configuration results in a weak attraction for the K⁺ ions in the bulk solution. The calculated local concentration of K⁺ (not shown) indicates that it is more probable to find K⁺ ions near the entrance and exit to the pore despite the fact that the channel carries a net positive charge. The repulsive influence of Arg⁸⁹ is counterbalanced by Asp⁸⁰ on the outer side, whereas the influence of Arg¹¹⁷ is counterbalanced by Glu¹¹⁸ on the inner side.

It is of interest to compare the calculated dynamical fluctuations with the information available from experimental data. In the absence of static disorder, the RMS fluctuations of the ith protein atom can be extracted from the crystallographic Debye-Waller B-factors (Brooks et al., 1988; Drenth, 1994; Willis and Pryor, 1975),

(1)

In Fig. 6, the RMS fluctuations calculated from the trajectory C1/MD-315K are compared with the x-ray data (R. MacKinnon, private communication). The statistical equivalence of the four subunits was exploited to improve the convergence of the calculated values. The quadratic fluctuations of atom i are written as a sum of two terms,

(2)

The first term, r_i²_(fluct), is the standard mean square atomic fluctuations averaged over the four monomers,

(3)

The second term, r_i²_(sym), represents the mean square deviation of the monomers away from the tetrameric symmetry,

(4)

where S_mr_i,m is the symmetrized average position of atom i in monomer m (S_m is a matrix operation, applied to monomers A, B, C, and D, corresponding, respectively, to rotations of 0, 90, 180, and 270° around the pore axis). The calculated fluctuations were averaged over the residues for consistency with the crystallographic structure refinement, which was performed with restrained B-factors (R. MacKinnon, private communication). A uniform offset constant of 0.6 Å was subtracted from the data to account for the presence of static disorder.

View larger version (18K): [in this window] [in a new window]   FIGURE 6   RMS fluctuations of the backbone atoms calculated from the trajectory at 315 K and from the experimental B-factors. All of the values are averaged over the individual amino acids. Equation 2 was used to exploit the tetrameric symmetry of KcsA in calculating the RMS fluctuations.

Although there are some important differences, the general form of the calculated fluctuations is in qualitative accord with the experimental data. The calculations show that the regions with the largest fluctuations are the NH₃⁺ and COO termini, the turrets (residues 50-60), and the short loop connecting the selectivity filter to the inner helix (residues 80-85). The main structural elements (outer, pore, and inner helices) have the smallest fluctuations.

The difference between the calculations and the data may be due to various factors. For example, it is likely that the fluctuations of the turrets are smaller in the x-ray data than in the calculations because this region of the channel is involved in protein-protein crystal contacts. On the other hand, the difference in the fluctuations of the protein chain near the NH₃⁺ terminus between the MD and the x-ray data could be due to the missing residues in the simulations. The simulated atomic system included only residues 23-119. However, the missing residues 1-22 and 120-126, which are present in the crystal even though they were not resolved in the crystallographic structure, can have an influence on the structure and dynamics of the rest of the protein. Last, alternative conformations of the unresolved residues at the NH₃⁺ terminus could give rise to inhomogeneous static disorder in the crystal.

Interestingly, a closer analysis reveals that the structure is systematically less flexible on the extracellular side than on the intracellular side. In Fig. 7, the RMS fluctuations are shown for all of the nonhydrogen atoms of KcsA as a function of their average Z position along the pore axis. A remarkable trend is observed. The value of the lowest RMS fluctuations is decreasing systematically in going from the intracellular side (Z < 0) to the extracellular side (Z > 0) of the structure. The trend involves both the backbone (black) as well as the side chain (red) atoms, although the backbone atoms usually have lower RMS. The selectivity filter is thus embedded in the most rigid part of the protein. The functional consequences of this dynamic feature of the channel are clear. The structure of KcsA is designed to have the less flexibility on the extracellular side on the membrane, where the channel must maintain the ability to discriminate between Na⁺ and K⁺ ions.

View larger version (32K): [in this window] [in a new window]   FIGURE 7   Calculated RMS fluctuations of all of the nonhydrogen atoms of KcsA as a function of their average Z position along the channel axis. The backbone (black) and the side-chain (red) atoms are shown.

Fluctuations of the selectivity filter

The selectivity filter, formed by the backbone carbonyl groups of residues Thr⁷⁵-Val⁷⁶-Gly⁷⁷-Tyr⁷⁸, is the key structural element responsible for the very high specificity of the K⁺ channel for K⁺ over Na⁺ ions. Its structure is both remarkable and simple: the main chain of each monomer runs almost vertically over a distance of 12 Å, with four carbonyl groups pointing in the same direction toward the center of a narrow pore (see Fig. 2). The structural stability of this functionally important region is significant. The calculated RMS fluctuations on the backbone atoms are on the order of 1.0 Å for this region of the channel. Examination of the crystallographic structure shows that a delicate lacework of hydrogen bonds plays an important role in maintaining the backbone of the residues forming the selectivity filter in a specific configuration. The most important hydrogen bonding pairs are given in Table 4 (see also Fig. 2).

At the bottom of the selectivity filter, the orientation of the carbonyl groups of Thr⁷⁴ and Thr⁷⁵ is stabilized by hydrogen bonds of the amide N-H group with the backbone carbonyl C==O of Thr⁷² and Glu⁷¹, respectively. Those hydrogen bonds appear to be stable on average, although they exhibit considerable fluctuations according to the current calculations.

Near the top of the selectivity filter on the extracellular side, the backbone amide group of Tyr⁷⁸ and Gly⁷⁹ and the carboxylate groups of Glu⁷¹ and Asp⁸⁰ are involved in a complex hydrogen-bonded network. The details of these hydrogen bonding interactions are somewhat uncertain because the side chain of Glu⁷¹ could not be completely resolved in the electron density, and its conformation is undetermined (Doyle et al., 1998). Its atomic coordinates are not included in the PDB file of the crystallographic structure. Furthermore, the ionization states of the carboxylate groups of Glu⁷¹ and Asp⁸⁰ are unknown. Previous MD simulations of the KcsA channel have assumed that those side chains are unprotonated (Guidoni et al., 1999; Shrivastava and Sansom, 1999). In the current calculations, Glu⁷¹ was modeled in a protonated (neutral) state to form a diacid hydrogen bond with the carboxylate group of Asp⁸⁰ near the extracellular surface, according to discussions with R. MacKinnon and co-workers (private communication). This choice is supported by a variety of evidence. In particular, a strong peak in the electron density map, observed near Asp⁸⁰, suggests that the side chain of Glu⁷¹ runs parallel to the pore axis toward the extracellular surface (R. MacKinnon, private communication). The peak of electron density, which remained unassigned in the refinement of the crystallographic structure, could correspond to the carboxylate group of Glu⁷¹; the peak can clearly be seen in the upper corner of figure 8 A of Doyle et al. (1998). A protonated state for Glu⁷¹ was chosen because the side chain is buried under the protein surface and is partly shielded from the bulk solvent. Its environment corresponds effectively to that of a low dielectric medium, and its pKa should be higher than normal. In support of the present decision, a valine at the corresponding position in Shaker (Val⁴³⁸) results in normal channel function, thus suggesting that substitution by a nonpolar side chain does not represent a significant disturbance on the structure. As shown in Table 4, the dicarboxylate Glu⁷¹ and Asp⁸⁰ form a strong hydrogen bond that remained very stable during all of the trajectories. The side chain of Glu⁷¹ is effectively shielded from the bulk solvent; there is only 0.5 water molecule on average within 3.4 Å of the protonated carboxylate group. Simultaneously, Asp⁸⁰ also forms a strong salt bridge with Arg⁸⁹ from a neighboring monomer. The salt bridge that was already present in the x-ray structure (although the distance between the charged side chains was initially larger; see Table 4) was rapidly strengthened and stabilized during all of the dynamic trajectories. A similar salt bridge was also observed in previous MD simulations of KcsA (Guidoni et al., 1999; Shrivastava and Sansom, 1999). A superposition of instantaneous conformations of the side chain of Glu⁷¹ is shown in Fig. 8. It is observed that alternative conformations occur, which may be one reason why the side chain could not be resolved in the electron density.

View larger version (14K): [in this window] [in a new window]   FIGURE 8   Superposition of alternative conformations of the side chain of Glu71 taken from simulation C1/MD-315K every 50 ps starting from t = 50 ps. The three populations of conformations correspond to different 1 angles separated by 120°.

Although the present model of Glu⁷¹ is reasonable, alternative solutions cannot be ruled out. For example, it is possible that the side chain is unprotonated and that the unassigned peak in the electron density corresponds to a water molecule (R. MacKinnon, private communication). Further calculations will be needed to fully resolve questions about the conformation and protonation state of Glu⁷¹.

The configuration of the selectivity filter near the bottom (Thr⁷⁴) and the top (Tyr⁷⁸) is stabilized by backbone amide N-H group hydrogen bonds. The situation is different for the carbonyl oxygen of Val⁷⁶, which is located in the center of the selectivity filter. Examination of the crystallographic structure shows that the orientation of the Val⁷⁶-Gly⁷⁷ peptide linkage is maintained only through a relatively weak hydrogen bond of the N-H group of Gly⁷⁷ with the backbone C==O group of Glu⁷¹ (see Fig. 2 and Table 4). The donor-acceptor distance is 3.25 Å in the x-ray structure, and the average calculated from the trajectories is around 3.0 Å. Interestingly, the residues Val⁷⁶ and Gly⁷⁷ have RMS fluctuations that are larger than those of their neighboring residues, according to both the simulation and the crystallographic data (see Fig. 6). Doyle et al. (1998) suggested that the larger fluctuations of Val⁷⁶ and Gly⁷⁷ might be due to multiple conformations of the backbone in response to alternative positions of the ions in the selectivity filter (see their note 19). Analysis of the simulations reveals that the amide plane formed by the peptide linkage Val⁷⁶-Gly⁷⁷ can undergo isomerization transitions, pointing alternatively into and away from the pore. The backbone dihedral (,)-map of the Val⁷⁶-Gly⁷⁷ amide plane is shown in Fig. 9 a. The backbone can adopt a configuration that corresponds alternatively to a right-handed -helix or a -sheet for Val⁷⁶, and a left-handed -helix or a -sheet for Gly⁷⁷. Snapshots of the two possible configurations are shown in Fig. 9 b. It is observed that the N-H group of Gly⁷⁷ can form a hydrogen bond with the water molecule located between the ions in the selectivity filter when the carbonyl of Val⁷⁶ points away from the pore. Because of the asymmetry of the water molecule, it is not possible to form favorable hydrogen bonds simultaneously with the carbonyl group of Val⁷⁶ of the four monomers. Further analysis shows that the transient isomerizations are also correlated with the movements of the K⁺ ions (see below).

View larger version (25K): [in this window] [in a new window]   FIGURE 9   (a) Phi-psi map of the Val76-Gly77 amide plane. Two configurations are illustrated. The backbone can adopt a configuration that corresponds alternately to a right-handed -helix ( = 60, = 60), and a -sheet ( = 90, = +120), for Val76, and a left-handed -helix ( = +60, = +60) and a -sheet ( = 90, = +120) for Gly77 (Schulz and Schirmer, 1979). (b) Snapshots of the two configurations of the carbonyl group of the Val76-Gly77 amide plane in the selectivity filter with the carbonyl pointing inward (left) or outward (right).

Several residues with aromatic side chains (Trp⁶⁷, Trp⁶⁸, Tyr⁷⁸, and Tyr⁸²) are disposed around the selectivity filter at the top of the KcsA structure on the extracellular side. They are well conserved in the sequence of other selective K⁺ channels, suggesting that they are functionally important. The packing of the aromatic side chains in KcsA is unusual. They lie parallel to the membrane plane, whereas aromatic side chains are generally packed in a perpendicular orientation in soluble proteins (Burley and Petsko, 1985). Doyle et al. (1998) suggested that the aromatic side chains contribute to stabilization of the channel structure around the selectivity filter, which is important for the high selectivity for K⁺ ions. To examine this hypothesis, the atomic fluctuations of the aromatic side chains in the directions parallel and perpendicular to the channel axis were calculated. An anisotropy factor  was defined for the atom i as the ratio of the fluctuations in the direction parallel to the channel axis, over the average fluctuations in the membrane plane,

(5)

The values were averaged over the four monomers of KcsA. The results, which are shown in Table 5, indicate that the atomic fluctuations are anisotropic. The structure of KcsA is thus significantly less flexible in the direction parallel to the membrane plane (i.e., perpendicular to the channel axis). Therefore, it seems likely that the role of the aromatic side chain of the residues Trp⁶⁷, Trp⁶⁸, Tyr⁷⁸, and Tyr⁸² is to strengthen the structure around the selectivity filter. This may be one of the factors responsible for the variations in the RMS fluctuations as a function of Z, which is observed in Fig. 7.

Dynamic transitions of ions

The high conductance of K⁺ channels is traditionally interpreted as an indication that they function as multiion pores (Hodgkin and Keynes, 1955; Hille and Schwarz, 1978). From a microscopic point of view, this implies that ion translocation is governed by concerted dynamic transitions involving several K⁺ ions simultaneously. It is thus of great interest to gain more information about such processes on the basis of atomic simulations. During the course of the first trajectory generated from configuration C1 at a temperature of 315 K (C1/MD-315K), the ions in the pore remained roughly near their initial position. To enhance the motion of the ions and water molecules inside the pore and increase the likelihood of observing dynamic transitions, a second, "hot" trajectory of 1.5 ns was generated at a temperature of 330 K. The initial configuration C1 was used again.

A concerted transition of the ions in the pore occurred after 250 ps during the C1/MD-330K trajectory. Observation of such rare events in a simulation, even though they may not be statistically significant, can give useful information. The time series of the Z coordinates of the ions and water molecules relative to the center of mass of the selectivity filter is shown in Fig. 10. The initial and final configurations of the ions and nearest water molecules in the selectivity filter are shown in Fig. 11. It is observed that the K⁺ ion that was initially in the cavity moves up into the narrow part of the selectivity filter at 250 ps, while the two other K⁺ ions move simultaneously toward the extracellular side. The final configuration, with three ions and two water molecules in the selectivity filter, remained stable for the rest of the trajectory. The entire transition took place over a period of ~400 ps. Although this is relatively short compared to physiological time scales, it is very long from the point of view of atomic motions in dense liquids. Thus the dynamic transition is not "ballistic" in character, but rather diffusive and overdamped, with many collisions between the ions, the water molecules, and the atoms forming the selectivity filter.

View larger version (56K): [in this window] [in a new window]   FIGURE 10   Concerted dynamical transition in the KcsA pore. The Z positions of the three K+ ions (green), a few selected water molecules (blue), and the oxygens of residues 75, 76, 77, and 78 (red) are shown as a function of time. All of the Z coordinates are given relative to the center of mass of the backbone of the residues forming the selectivity filter. The positions of the carbonyls were averaged over the four subunits. The transition occurs over a period of ~400 ps, from t = 250 ps to t = 650 ps. At t = 250 ps the Val76-Gly77 amide plane of one subunit undergoes an isomerization transition from the outward to an inward configuration (see Fig. 9).

View larger version (36K): [in this window] [in a new window]   FIGURE 11   Initial (left) and final (right) states of the ions and water molecules in the pore from simulation C1/MD-330K at t = 0 ps and t = 1.5 ns, respectively.

The most stable configurations (i.e., before and after the transition) correspond to a single file of ions separated by a single water molecule, with two or three K⁺ occupying the selectivity filter, with one (K-W-K) or two (K-W-K-W-K) water molecules, respectively. A water molecule is also strongly bound to the ion located at the inner site on the cavity side. On a speculative note, it is conceivable that the configurations with two and three K⁺ ions correspond to the discrete states postulated in kinetic rate models of multiion channels (Hodgkin and Keynes, 1955; Hille and Schwarz, 1978).

The number of nearest neighbors (water and carbonyl oxygens), shown in Fig. 12, varies in a complex nonmonotonic way during the dynamic transition. Before the transition, the K⁺ ion in the outer site is coordinated by the eight carbonyl oxygens of Gly⁷⁷ and Tyr⁷⁸ of the four monomers and makes contact with no water molecules, the K⁺ ion in the inner site is coordinated by the carbonyl oxygens of Thr⁷⁵ and Val⁷⁶ and makes contact with two water molecules, and the K⁺ ion in the cavity is solvated by six to eight water molecules. After the transition, the coordination shell of the K⁺ ions is significantly transformed. All three ions are now in contact on average with two water molecules. The ion in the outer site has moved up along the channel axis and is now surrounded by the carbonyl oxygen of Tyr⁷⁸, while the ion in the inner site is surrounded by the carbonyl oxygens of Val⁷⁶ and Gly⁷⁷. The ion that was in the cavity has entered the selectivity filter and is surrounded by the carbonyl oxygens of Thr⁷⁵. This ion also makes some transient contacts with the OH group of the side chain of Thr⁷⁵ of one of the monomers.

View larger version (46K): [in this window] [in a new window]   FIGURE 12   Variation of the number of water and carbonyl nearest neighbors during the transition for each of the ions. Throughout the trajectory all three ions are coordinated by a total of six to eight water molecules or carbonyls. After the transition, the lowest ion also makes some contacts with the OH group of the side chain of Thr75 of one of the monomers (not shown).

The dynamic fluctuations of the channel structure and the movements of the ions appear to be coupled. The orientation of the Val⁷⁶-Gly⁷⁷ peptide linkage is affected by the position of the ions in the selectivity filter. When there are two ions in the selectivity filter, it is observed that the carbonyl group of Val⁷⁶-Gly⁷⁷ of one of the four subunits of the KcsA tetramer undergoes transient isomerization transitions, sometimes pointing away from the center of the pore (see Fig. 9), while no such isomerizations are observed once the selectivity filter is occupied by three ions. Last, it is observed that several microscopic events are correlated with the dynamic transition of the ions. Those involve the reorientation of the carbonyl group of Val⁷⁶ (see Fig. 9 and caption of Fig. 10) as well as the fluctuations in the hydrogen bond of Gly⁷⁷ N-H with the backbone C==O of Glu⁷¹, the dicarboxylate hydrogen bond between the side chains of Glu⁷¹ and Asp⁸⁰, the hydrogen bond between the aromatic residues Trp⁶⁸ and Tyr⁷⁸, and the salt bridge between Asp⁸⁰ and Arg⁸⁹ (see Fig. 2). The distance between these molecular moieties is on the order of 10-15 Å. Therefore, ion translocation in the pore appears to be coupled to structural fluctuations of the channel occurring at very large distances from the ions.

K⁺ ions and water molecules in the selectivity filter

So far, we have examined the trajectories generated from the initial configuration C1 with one K⁺ ion in the outer site, one K⁺ ion in the upper-inner site, and a single water molecule in the selectivity filter (see Table 1). This initial configuration is based on the atomic coordinates of the ions and water molecule of the crystallographic structure. Although this configuration is certainly expected to play an important role in the conduction of K⁺ ions through the KcsA channel, it is also of interest to examine other configurations of the ions and water molecules in the selectivity filter.

Ion conduction is intrinsically a dynamic process in which a number of different configurations may play an important role. In addition, the possibility of alternative configurations of the ions and water molecules in the selectivity filter, which would be consistent with the available experimental data, has been suggested by Doyle et al. (1998). The positions of the K⁺ ions reported in the crystallographic structure were determined experimentally by analyzing difference Fourier maps of crystals containing Rb⁺ and Cs⁺ relative to those with K⁺. Excess positive electron density revealed the position of K⁺ ions replaced by the more electron-dense Rb⁺ or Cs⁺ ions in the outer and inner sites. Two peaks of excess density observed in the difference Fourier map obtained at 4.0-Å resolution with Rb⁺ ions were interpreted as alternative positions for the inner site (the "upper" and "lower" inner site; see Fig. 2). Diffuse electronic density was also detected in the center of the cavity, although no ion was explicitly included at this position in the crystallographic structure. Difference Fourier analysis of ion-substituted crystals indicated that a hydrated monovalent cation is present in the cavity (Roux and MacKinnon, 1999). A strong peak of electron density located between the inner and outer sites in the selectivity filter did not vary in the difference Fourier maps and was attributed to a single water molecule. Nonetheless, the exact number of water molecules between the ions in the selectivity filter is uncertain. The distance between the ions in the outer and the upper inner sites is 7.3 Å, and the water molecule is at a distance of 3.4 Å and 3.9 Å from the ions. Such distances appear to be slightly large compared to the optimal distance between a K⁺ and a water molecule. Such a distance is ~2.6 Å according to ab initio calculations (Roux and Karplus, 1995). Although a single water molecule was included in the crystallographic structure, additional water molecules located between the ions in the selectivity filter could have remained undetected experimentally at the current resolution. Accurate determination of solvent molecule positions in a protein crystallographic structure often requires a resolution on the order of 2.0 Å (Drenth, 1994; Burling et al., 1996).

Three alternative configurations of the ions and water molecules in the selectivity filter were considered. The three configurations (C2, C3, and C4) were constructed from the atomic coordinates taken from C1 after 300 ps of equilibration (see Fig. 4). Those configurations were then refined using the following procedure: the ions in the narrow part of the selectivity filter were held fixed in place using strong harmonic energy restraints, and all of the atoms within a distance of 6 Å were relaxed and optimized using energy minimization and MD. Details of the ions and water molecules in all of the initial configurations are given in Table 1. Each of the configurations is consistent with the available experimental data (Doyle et al., 1998) and is possible on structural grounds. For example, in C2 the ion-water distances are 2.6 Å and the water-water distance is 2.7 Å. In C4 the ion-water distances are 2.6 and 2.7 Å and the water-water distances are around 2.8 to 2.9 Å. Nonetheless, further refinement of the three additional configurations shows that only C3 is stable. During the structural refinement of C2 and C4, one of the water molecules in the selectivity filter was rapidly expelled through the wall of the channel. Although the extra water molecule does not make any obvious bad contact with another atom, it is unstable because it cannot achieve favorable electrostatic interactions with the two ions in the narrow filter. The unrealistic pathway followed by the extra water molecule, which is probably caused by the excessive repulsive forces, clearly indicate that configurations C2 and C4 are not possible. In contrast, the configuration C3 (outer/lower-inner with two waters) could be relaxed without large distortions of the channel and displacements of ions and water molecules. To further examine its stability, two trajectories of 0.5 ns at 315 K (C3/MD-315K) and 330 K (C3/MD-330K) were generated following the same procedure as that used for the trajectories C1/MD-315K and C1/MD-330K. The C3 configuration converted rapidly into a C1-like configuration in less than 300 ps during the equilibration with progressively decreasing restraints. Fig. 13 shows the initial and final states of the ions and water molecules in the selectivity filter during trajectory C3/MD-330K. This computer experiment, therefore, shows that the configuration C3, with two ions separated by two water molecules, is not a stable state of the system.

View larger version (24K): [in this window] [in a new window]   FIGURE 13   Initial (left) and final (right) states of the trajectory C3/MD-300K, which was generated from the configuration C3 (see Table 1). The initial configuration was taken at the beginning of the equilibration of configuration C3; the final one was taken from simulation C3/MD-330K at t = 500 ps.

On the basis of the trajectories C1/MD-330K, C3/MD-315K, and C3/MD-330K, we conclude that there are only two stable states for the K⁺ ions and water molecules in the selectivity filter (even though their relative free energies and populations cannot be determined by the current calculations). A first state corresponds to two ions (K) located in the selectivity filter separated by a single water (W) molecule, K-W-K. A second state corresponds to three ions in the selectivity filter with one water molecule separating the ions, K-W-K-W-K. In both states, a single water molecule is stabilized between the two ions of each pair, suggesting that the K⁺ cation translocating though the narrow pore is accompanied by exactly one water molecule. This observation is in qualitative accord with streaming potential measurements, which suggest that one or two water molecules per ion flow through the channel (Alcayaga et al., 1989).

In fact, the two stable states are structurally related to the cation binding sites and crystallographic water molecule that were detected experimentally. As shown schematically in Fig. 14, there are five stable positions or "sites" in the selectivity filter that can be occupied by a water molecule or a K⁺ ion. The first (K-W-K) and second (K-W-K-W-K) states are observed in the course of the trajectories C1/MD-330K, C3/MD-315K, and C3/MD-330K. On a speculative note, it is possible that the two states with two and three ions are superimposed in the KcsA crystals with K⁺ ions. One implication of this hypothesis is that the positions of the maxima observed in the difference Fourier maps result from the combined influence of slight variations in the position of the substituted ions and differences in the populations of the two stable states relative to K⁺. One may expect slight variations in the relative free energies of the two- and three-ion states with Rb⁺, K⁺, and Cs⁺ because of ion-ion interactions (Roux et al., 1995). Further computational studies of the KcsA channel with Rb⁺ and Cs⁺ may lead to better interpretation of the experimental data for substituted ions.

View larger version (23K): [in this window] [in a new window]   FIGURE 14   Schematic representation of the two stable states of the K+ and water molecules in the selectivity filter. There are five local sites. In the first state (left), the two K+ ions occupy the crystallographic upper inner (4) and outer (2) site, respectively. One water molecule occupies the position assigned to the crystallographic K+ ion in the lower inner site (5), and one water molecule occupies the position assigned to the crystallographic water molecule (3). No water molecule is located above the K+ in the outer site (i.e., site 1 is empty in this state). In the second stable state (right), one K+ ion occupies the crystallographic lower inner outer site (5), one K+ ion occupies the position assigned to the crystallographic water molecule (3), and one K+ ion is located above the crystallographic outer site (1). One water molecule occupies the position assigned to the crystallographic K+ ion upper inner site (4), and one water molecule occupies the position assigned to the crystallographic K+ ion outer site (2).

The observed configuration with three K⁺ and water molecules in single file seems to challenge our view of ion-ion repulsion. Can three K⁺ ions be so close in space without repelling one another? Or perhaps, one may more soberly ask whether a configuration with three K⁺ in the selectivity filter is plausible at all or is an artifact of the calculations. Clearly, further calculations will be required to quantitatively characterize the configurations of the K⁺ ions inside the KcsA channel. Nonetheless, in support of the current observations one may note that the TIP3P water potential (Jorgensen et al., 1983)), which does not account for induced electronic polarization, was used for simulations. Because of the lack of induced polarization, the ion-ion repulsion is undershielded in the present simulations. Therefore, the configuration with three K⁺ ions in the selectivity filter that was observed, despite an underestimated shielding of the ion-ion repulsion, might be even more stabilized if a potential function allowing for induced polarization had been used for the simulations.

	CONCLUSION

TOP ABSTRACT INTRODUCTION THEORY AND METHODS RESULTS AND DISCUSSION CONCLUSION REFERENCES

In this paper we described the results from MD simulations of an atomic model of the KcsA K⁺ channel embedded in an explicit membrane bathed by a 150 mM KCl aqueous salt solution. In all of the simulations, the three-dimensional structure of KcsA in the membrane system was very stable. In particular, the overall RMS deviation of the whole tetramer is on the order of 1.9 Å relative to the crystallographic structure, while the deviation of all of the main secondary elements is less than 0.9 Å. Such a stability is a strong indication that the fold and packing of KcsA are accurate despite the moderate resolution of the x-ray data.

Although the current calculations show that the atomic fluctuations are relatively small in the functionally important regions of the structure, their magnitude is significantly larger than the difference in radius between a Na⁺ ion and a K⁺ ion. The RMS fluctuations are on the order of 1.0 Å for the residues forming the selectivity filter. The hydrogen bonds, which play an important role in maintaining the selectivity filter in its configuration, are stable on average, although they exhibit considerable fluctuations. Residues with aromatic side chains (Trp⁶⁷, Trp⁶⁸, Tyr⁷⁸, and Tyr⁸²) help to strengthen and stabilize the channel structure around the selectivity filter. Their atomic fluctuations are anisotropic, indicating that the structure of KcsA is significantly less flexible in the direction perpendicular to the channel axis. The fluctuations of Val⁷⁶ and Gly⁷⁷ are slightly larger, in agreement with the experimental B-factors, because the amide plane can undergo isomerization transit