INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS REFERENCES

Ion channels play a key role in the electrical activity of excitable cells (Hille, 1992), enabling passive movement of ions across membranes. Ions move at high rates (~10⁷ ions s¹ channel¹), and yet channels can be highly selective as to which ions may pass. Furthermore, channels are gated, i.e., they open and close in response to changes in transmembrane voltage and/or ligand binding to the channel protein. In the membranes of excitable cells potassium (K) selective channels are responsible for the repolarizing phase of action potentials and for controlling membrane excitability. Furthermore, potassium channels play diverse roles in a wide range of cells, in organisms ranging from bacteria to plants and animals. In particular, certain potassium channels (so-called "background" K channels; Maingret et al., 1999) seem to be responsible for the passive permeability of cell membranes to K⁺ ions, which is central to generation of a voltage difference across cell membranes. Despite this diversity of function, potassium channels seem to share a common pore-lining domain, composed of four repeats of a motif made up of two transmembrane (TM) helices flanking a re-entrant P-loop which carries the main determinants of ion selectivity (MacKinnon et al., 1998; Miller, 1991).

The structure of a bacterial potassium channel (KcsA from Streptomyces lividans; Schrempf et al., 1995) has been solved by x-ray diffraction at 0.32 nm resolution (Doyle et al., 1998). Given the conservation of the pore-lining domain between the various potassium channels, this structure provides a framework within which to understand K⁺ selectivity and permeation. However, a crystal structure inevitably provides a static, spatially and temporally averaged image of a channel. To bridge the gap between molecular structure and physiological behavior an understanding of the atomic resolution dynamics of potassium channels is required. One way in which to approach this is via simulation studies. This approach complements, continuum electrostatics calculations (Roux and MacKinnon, 1999) by providing a more dynamic image of channel function.

Molecular dynamics (MD) simulations enable one to explore the motions of biomolecules on a picosecond-to-multinanosecond time scale (Brooks et al., 1988; Daura et al., 1998; Duan and Kollman, 1998; McCammon and Harvey, 1987; van Gunsteren and Mark, 1992). In particular, it is now possible to simulate fully solvated lipid bilayers in a realistic fashion (Jakobsson, 1997; Tieleman et al., 1997; Tobias et al., 1997). Simulations of pure lipid bilayers have been extended to simulations of lipid bilayers containing TM peptides (Belohorcova et al., 1997; Bernèche et al., 1998; Biggin and Sansom, 1998; Chiu et al., 1999; Forrest et al., 1999; Roux and Woolf, 1996; Shen et al., 1997; Tieleman et al., 1999c; Woolf and Roux, 1994), assemblies of TM peptides (Tieleman et al., 1999a), or large integral membrane proteins (Tieleman and Berendsen, 1998). These simulations allow one to explore the conformational dynamics of membrane-spanning peptides and proteins, and their interactions with their environment (Tieleman et al., 1999b) in some detail.

MD simulations of channel proteins also provide considerable information on channel/ion/water interactions. Pioneering studies that applied this approach to model channels formed by gramicidin (Chiu et al., 1996, 1999; Roux and Karplus, 1994) have since been extended to channels formed by -helix bundles such as alamethicin (Tieleman et al., 1999a) or synthetic channel-forming peptides (Randa et al., 1999; Zhong et al., 1998b, c). To fully exploit the power of simulations to understand the energetics and dynamics of ion channels it is important to include as complete a representation as possible of the anisotropic environment provided by a lipid bilayer. A fist approximation to simulation in a lipid bilayer may be obtained via simulation in an octane "slab," which is solvated on either side and into which the membrane protein is inserted. This approach has been used by Klein and colleagues to study a number of ion channels formed by bundles of peptide helices (Zhong et al., 1998a-c) and more recently has been extended to KcsA (Guidoni et al., 1999). However, it is feasible to represent a lipid bilayer explicitly in such simulations (Chiu et al., 1999). This is important, as both the presence of lipid headgroups and the fluidity properties of a lipid bilayer (which may differ significantly from those of an isotropic solvent such as octane) may influence the dynamic properties of the embedded channel. Indeed, some such differences have been suggested by comparison of the behavior of simple peptide channels simulated in an octane slab (Zhong et al., 1998b, c) and in a phospholipid bilayer (Randa et al., 1999).

In this paper we present several simulations of a potassium channel (KcsA) embedded in a phospholipid bilayer. In our analysis of these simulation we focus on movements of K⁺ ions and water through the central pore. Each of our five simulations is of 1 ns duration. This is about an order of magnitude shorter than the mean time for passage of an ion through a channel (~15 ns for a 100 pS conductance channel at 100 mV transbilayer voltage). Thus one may hope to capture some aspects of ion permeation, although longer simulations will be required for proper statistical sampling. Indeed, to fully address all aspects of ion permeation it may be necessary, in addition to MD simulations, to run longer-time scale coarse-grained simulations using methods such as Brownian dynamics (Bek and Jakobsson, 1994).

The structure of KcsA is that of a truncated cone, with a central pore running down the center (Fig. 1). The wider end of the cone corresponds to the extracellular mouth of the channel. This contains a selectivity filter that is formed by a TVGYG sequence motif characteristic of potassium channels. In the x-ray structure this adopts an irregular conformation, and binds two K⁺ ions plus an intervening water molecule. Beneath the selectivity filter is a central water-filled cavity, which also appears to contain a (loosely bound) cation in the x-ray structure. Finally, the pore constricts at its intracellular mouth to form a putative gate region where the pore radius falls to ~0.11 nm (i.e., less than the Pauling radius of a K⁺ ion, 0.13 nm). Our simulations reveal a concerted movement of water molecules and K⁺ ions within the selectivity filter, and suggest how breathing motions of the channel protein (on a ~0.1-ns time scale) may transiently open the intracellular gate, allowing a K⁺ ion to leave the channel.

View larger version (48K): [in this window] [in a new window]   FIGURE 1   Overview of simulations. (A) Location of the ion sites within the channel structure. On the left-hand side a diagram of the selectivity filter (T75VGY78) is given (for simplicity, only two subunits are shown) along with the ion (and water) binding sites defined in the x-ray structure. Sites S1-S4 form the selectivity filter region and are defined by the O atoms. The cavity site C is also shown. The right-hand side of the diagram defines the locations of the K+ ions and crystallographic water molecule included in the initial configurations of the simulations. Note that those sites not occupied by a K+ ion or crystallographic water may become occupied by water molecules during the solvation stage of the simulation setup procedure. The ions are named K1, K2, and K3 according to their initial locations. This nomenclature will be used throughout the discussion of the simulations. (B) Three snapshots of simulation MDK3. Two subunits (B and D) of the four that make up the KcsA channel are shown, plus the three K+ ions (K1, K2, and K3). The M1 helix, P-region, and M2 helix are shown as red, green, and blue ribbons, respectively. The backbone of the selectivity filter and the side chains of T75 and T107 are shown as individual bonds between atoms. The horizontal white dotted lines indicate the approximate position of the phospholipid headgroups, and the vertical blocks to the left of the diagram indicate the extent (along the pore axis z) of the selectivity filter (s; green), central cavity (c; gray), and gate (g; blue) regions. The snapshots correspond to the start of the simulation (100 ps), late in the simulation just before ion K3 leaves the pore (850 ps) and at the end of the simulation when K3 is on the intracellular side of the membrane (1100 ps). Comparing 100 ps and 850 ps one can see the concerted movement of K1 and K2 between sites of the selectivity filter, and the movement of K3 deeper into the cavity to interact with one of the T107 side chains. Comparing 850 ps and 1100 ps one can see the gradual movement of K2 toward the cavity and the exit of K3 through the intracellular mouth of the pore.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS REFERENCES

Simulation systems

All simulations started from the same protein structure and bilayer model. The protein structure was that in PDB file 1bl8, with the modification that atoms for all side chains were included. The "missing" residues were added by building in stereochemically preferred conformers. Although not present in the x-ray structure, an inter-subunit salt bridge (D80 to R89) was formed during the early stages of the simulations. Note that residues 1-22 and 120-160 are missing from this model, which therefore represents the core channel-forming domain of the protein. A fully equilibrated palmitoyl oleoyl phosphatidylcholine (POPC) lipid bilayer (kindly provided by Dr. Peter Tieleman) was used as the starting point for generating the phospholipid bilayer into which KcsA was embedded (see below).

Five simulations were performed, differing in the number and locations of K⁺ ions in the initial model. In describing these models we will use the ion/water site nomenclature defined in Fig. 1 A. Thus, sites S1 to S4 make up the selectivity filter and site C is that near the center of the cavity. In the x-ray structure, sites S1 and S3 or S4 are occupied by a cation, while S2 is occupied by a water molecule. In the subsequent discussions this water will be referred to as the "crystallographic" water molecule, W1. In simulation MDK0 no K⁺ ions were included, either in the selectivity filter or in the cavity. In MDK1, a single K⁺ ion (henceforth referred to as K1) was included at site S1 (Fig. 1). In MDK2, two K⁺ ions (K1 and K2) were included at sites S1 and S3, with the crystallographic water at S2. In simulation MDK2' the second ion (K2) occupied site S4. Thus MDK2 and MDK2' correspond to the two alternative configurations of ion found in the x-ray structure, in which sites S3 and S4 exhibit partial occupancy (Doyle et al., 1998). Simulation MDK3 started from the same structure as MDK2, except for a further ion (K3) in the cavity.

Setup of simulations

The simulation methodology was similar to that used in MD simulations of the bacterial porin OmpF (Tieleman and Berendsen, 1998), of a channel formed by alamethicin (Tieleman et al., 1999a), and of channels formed by a synthetic peptide (Randa et al., 1999), all of which have yielded reasonable correlations with experimental data. In the latter two studies a pre-existing equilibrated POPC bilayer was used, into which a cylindrical hole was introduced, by a combination of removing a small number of lipid molecules and running a short MD simulation with a radially acting repulsive force in order to drive any remaining lipid atoms out of the cylinder into the bilayer. The channel molecule was then embedded within the hole thus created. A similar protocol was adopted in this study. However, the more asymmetric shape of KcsA required that after embedding the protein, a somewhat longer preparatory MD simulation was needed to re-pack lipid molecules around the protein (see below). Thus, the KcsA channel (containing four subunits each of 97 residues) was embedded in a lipid bilayer containing 243 POPC molecules (116 in the upper, i.e., extracellular, leaflet and 127 in the lower, i.e., intracellular, leaflet). This provided a difference in lipid surface area between the two leaflets approximately equivalent to the difference in cross-sectional area of KcsA at the extracellular and intracellular end of the molecule. The ions and/or crystallographic water were then added. The resultant system was solvated with ~10,000 SPC water molecules, giving a total of ~42,000 atoms. Upon solvation we observed the cavity to contain from 15 to 18 water molecules. This was in agreement with a calculation of the cavity volume based on integration of the pore radius profile, and was robust to either small changes in van der Waals radii of the waters or to use of either standard GROMACS or CHARMM solvation procedures (Shrivastava and Ranatunga, unpublished results). All ionizable residues were in their default protonation state. Sufficient Cl ions (from four (MDK0) to seven (MDK3)) were added to the bulk solvent on either side of the bilayer to give a net charge of zero. These ions replaced water molecules at the positions of lowest Coulomb potential. This was done by removal of successive water molecules, one at a time, and calculation of the Coulombic interaction energy of a Cl ion at that position with the remainder of the system.

Following setup of the system close packing of the protein molecule and the phospholipids was achieved by a 100-ps simulation, during which the protein atoms were fixed and a lateral pressure of 500 bar was applied. Having thus embedded the protein in the bilayer, the protein and potassium ion coordinates were restrained during a further 100-ps equilibration period during which no excess lateral pressure was applied and lipids and water were free to move. Finally, all restraints were removed during the subsequent 1-ns production runs.

Simulation methodology

MD simulations were run using GROMACS (http://rugmd0.chem.rug.nl/~gmx/gmx.html). A twin range cutoff was used for longer-range interactions: 1.0 nm for van der Waals interactions and 1.7 nm for electrostatic interactions. The time step was 2 fs, with the LINCS algorithm to constrain bond lengths. We used NPT conditions in the simulation. A constant pressure of 1 bar independently in all three directions was used, with a coupling constant of _p = 1.0 ps (Berendsen et al., 1984). Water, lipid, and protein were coupled separately to a temperature bath at 300 K, using a coupling constant _T = 0.1 ps. MD simulations were performed on an 80-node SGI Origin 2000, typically taking ~10 days of cpu time on eight R10000 processors.

The lipid parameters were as in previous MD studies of DPPC bilayers (Berger et al., 1997; Marrink et al., 1998) (with the addition of some GROMOS parameters for the double bond in the acyl tail), and as in previous MD simulations of Alm (Tieleman et al., 1999a, c). These lipid parameters give good reproduction of the experimental properties of a DPPC bilayer. The lipid-protein interactions used GROMOS parameters. The water model used was SPC (Berendsen et al., 1981), which has been shown to be a reasonable choice for lipid bilayer simulations (Tieleman and Berendsen, 1996). The K⁺ parameters (kindly supplied by Dr. Peter Tieleman) were as in Straatsma and Berendsen (1988).

Structural diagrams were prepared using Molscript (Kraulis, 1991) and Raster3D (Merritt and Bacon, 1997). Pore radius profiles were calculated using HOLE (Smart et al., 1993). Secondary structure analysis used DSSP (Kabsch and Sander, 1983). Other analyses used GROMACS and/or locally written code.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS REFERENCES

Drift, fluctuations, and secondary structure

Analysis of the all-atom RMSD (root-mean-square deviation) versus time revealed a similar pattern for each of the five simulations, namely an initial jump (over ~100 ps) to an RMSD of 0.2 nm, followed by little further drift over the next few hundred picoseconds to a plateau value of ~0.25 nm, which is maintained throughout the rest of the simulation (Fig. 2). This is typical of MD simulations of a membrane protein in a bilayer (see, e.g., Tieleman and Berendsen, 1998). The initial jump in RMSD is presumed to reflect relaxation of the protein upon transfer from a crystal to a bilayer environment and/or inaccuracies in the potential function. There was no significant difference in overall RMSD among the five different simulations, suggesting that the presence/absence of K⁺ ions does not influence the overall conformational stability of the protein.

View larger version (32K): [in this window] [in a new window]   FIGURE 2   Drift of protein structure from the initial model. The root-mean-square deviation (RMSD) of all protein atoms of KcsA from the starting structure is shown as a function of time for all five simulations. (Note that these and subsequent graphs start at 100 ps, as during 0 to 100 ps the protein coordinates were restrained.)

Identification of the more flexible regions of a protein during a simulation may be obtained via examination of the root-mean-square fluctuation (RMSF) of the C atom of each residue from its time-averaged position (Fig. 3). Similar overall patterns of RMSF versus residue number were seen for all five simulations, with the fluctuations ranging from 0.05 to 0.25 nm. The peaks in the C RMSF values are generally observed at the N- and C-termini of each subunit and in the extracellular surface loops on either side of the P-region. These residues have the largest B-factors in the crystal structure. In particular, the loop between the M1 helix and the P helix (i.e., the "turret" loop) seems to fluctuate most markedly. As this loop is oriented away from the rest of the protein, toward the surrounding solvent, it is not unreasonable that it shows the greatest mobility. Comparison of the five simulations shows no major difference. All exhibit C RMSF values of ~0.05 nm for the cores of the M1 and M2 helices. However, there is some suggestion of greater fluctuations in the selectivity filter residues for simulations MDK0 and MDK1 compared to MDK2, MDK2', and MDK3 (see below).

View larger version (40K): [in this window] [in a new window]   FIGURE 3   Fluctuations of protein coordinates. The root-mean-square fluctuations (RMSFs) of the C coordinates from their time-averaged values are shown as a function of residue number for all five simulations. The dotted vertical lines indicate the divisions among the four subunits. Note that the residue numbers are relative, such that residue 1 corresponds to residue 22 of subunit A in the x-ray structure, residue 98 corresponds to residue 22 of subunit B, residue 195 corresponds to residue 22 of subunit C, and residue 292 corresponds to residue 22 of subunit D. Thus the selectivity filter regions (indicated by small vertical arrows) correspond to residues 53-56, 150-153, 247-250, and 344-347 in this and subsequent figures.

The consequences of these fluctuations on the secondary structure of KcsA can be visualized using DSSP (Kabsch and Sander, 1983). As can be seen from Fig. 4, the overall secondary structure pattern of KcsA is maintained, but changes occur in the loop regions as a function of time. The selectivity filter region in MDK1, MDK2, MDK2', and MDK3 is seen to be relatively stable, whereas some transient changes in secondary structure are seen for the same region in MDK0. In the case of MDK3, there is limited loss of -helicity in the center of M2 (in the vicinity of residues 174-177), in subunit B. When visualized, this reveals a slight uncoiling of the M2_B helix, accompanied with a slight inward bend.

View larger version (46K): [in this window] [in a new window]   FIGURE 4   Secondary structure (analyzed using DSSP; Kabsch and Sander, 1983) of KcsA as a function of time for simulations MDK0, MDK1, MDK2, and MDK3. Black indicates -helical residues, dark gray 310-helix, light gray turn, and white indicates random coil. Note that residue numbering is the same as in Fig. 3. The horizontal arrow indicates the distorted region of helix M2B in simulation MDK3. The M1, P-helix, filter (f), and M2 helix are indicated alongside the plot for MDK0.

K⁺ ions stabilize the selectivity filter

The selectivity filter is composed of five residues, T⁷⁵VGY⁷⁸G. As it has been suggested, on the basis of physiological studies of potassium channels, that K⁺ ions may be required to maintain the structural integrity of the selectivity filter (Khodakhah et al., 1998; Melishchuk et al., 1998; Ogielska and Aldrich, 1999), we have examined the structural dynamics of this region in more detail. The RMSD values of all atoms in this region are compared for the five simulations in Fig. 5. In MDK0 and MDK1, the RMSD is from 0.17 to 0.19 nm. In contrast, in simulations MDK2, MDK2', and MDK3 it is significantly lower (~0.11 nm in MDK3 until the last 100 ps of the simulation). This suggests that the simultaneous presence of two K⁺ ions within the selectivity filter tends to stabilize its conformation, whereas in the presence of a single ion or no ions, conformational changes may occur. Significantly, there is a jump in RMSD for the selectivity filter in simulation MDK3 at ~1000 ps, corresponding to the exit of K2 from the filter (see below) leaving K1 at site S2 as the only ion in the filter.

View larger version (34K): [in this window] [in a new window]   FIGURE 5   Drift of the selectivity filter structure from the initial model. The RMSD of all atoms of the selectivity filter (T75VGY78) from its starting structure is shown as a function of time for all five simulations.

To explore these conformational changes in more detail we examined Ramachandran plots of the selectivity filter residues as a function of time. These confirm that there is significant flexibility in the selectivity filter backbone in the absence of K⁺ ions (simulation MDK0), which is largely suppressed in the presence of K⁺ ions (e.g., simulation MDK3). For example, in Fig. 6 the Ramachandran plots of residues G77 and Y78 of a single subunit from MDK0 and MDK3 are compared. In MDK0 both G77 and Y78 change conformation such that, after ~300 ps, they differ from the x-ray structure. In contrast, in MDK3 both residues exhibit limited fluctuations about the initial conformation. This flexibility in MDK0 was observed in the other three subunits, and also in simulation MDK1. MDK2' exhibited much the same behavior as MDK2. Thus, it seems that at least two K⁺ ions need to be present within the selectivity filter to maintain the conformation seen in the crystal structure.

View larger version (21K): [in this window] [in a new window]   FIGURE 6   Backbone torsion angle dynamics in the selectivity filter. Ramachandran plots of residues G77 and Y78 in simulations MDK0 and MDK3. In each case the plot corresponds to a single subunit and shows the torsion angles at all times during the simulation. The arrows in the MDK0 plot indicate the change in torsion angles as time progresses.

In the remainder of the paper we will therefore focus our analysis on those simulations (i.e., MDK2, MDK2', and MDK3) in which the selectivity filter remains intact. In particular we will concentrate on simulation MDK3, as this corresponds most closely, in its initial configuration, to the x-ray structure. Furthermore, at ~900 ps ion K3 leaves the channel cavity, via the intracellular mouth (Fig. 1 B), and ion K2 starts to move from the selectivity filter into the cavity. Thus MDK3 provides information on the dynamic events underlying ion permeation through KcsA.

Concerted movement of ions and water in the filter

Movements of the K⁺ ions and crystallographic water molecule in the selectivity filter were examined in terms of their z-coordinates as a function of time (Fig. 7). At first site it appears that K1, W1, and K2 simply remain in the filter throughout all three simulations. Closer examination reveals a subtler pattern, corresponding to a concerted, single-file motion of ions and water. For example, in MDK3 at ~350 ps ions K1 (initially at site S1) and K2 (initially at site S3) both move one site down the selectivity filter to sites S2 and S4, respectively. These movements are in concert with one another, and with that of the intervening crystallographic water from site S2 to S3. In simulation MDK2, the movements of the two ions are similar to those of K1 and K2 in MDK3, but occur at ~150 ps. In MDK2' there is a comparable transition at ~830 ps, when K1 moves from S1 to S2 and the water from S2 to S3. (Note that before this transition there is a "vacuum" between the water and K2. It is possible that in the x-ray structure, when an ion is at S4, there is a second water at S3. This is being explored in further simulations; Shrivastava and Sansom, unpublished results). In contrast, in MDK1 (data not shown), the single K⁺ ion makes an initial jump from site S1 to S2, where it then remains for nearly 1 ns. Snapshots of MDK3 at early and late stages in the simulation (Fig. 8) show that in addition to the crystallographic water initially at site S2, a water molecule initially at the extracellular entrance to the selectivity filter moves to site S1 to replace K1. Thus, there is concerted movement of a water-K⁺-water-K⁺ column through the filter on a time scale of several hundred picoseconds.

View larger version (43K): [in this window] [in a new window]   FIGURE 7   Ion/water z trajectories for simulations MDK2, MDK2', and MDK3. The time-dependent locations on z of K1 and K2 are shown, along with that of the crystallographic water (W1, gray lines) that is initially located at site S2. The vertical arrows show the approximate times of the concerted K1-water-K2 movement within the filter.

View larger version (31K): [in this window] [in a new window]   FIGURE 8   The selectivity filter at 100 ps and 850 ps in simulation MDK3, showing the concerted motion of ions and water. At 100 ps there is a water molecule (dark gray/white) just extracellular to the filter (top of diagram), K+ ions (pale gray) at sites S1 and S3, and water molecules at site S2 and just below site S4. At 850 ps there are water molecules at sites S1 and S3 and just below site S4, and K+ ions at sites S2 and S4. The peptide backbone of the filter and the T75 side chains are shown (for subunits A and C only, for clarity).

Examination of Fig. 8 also suggests subtle changes in backbone conformation to maintain optimal interactions between the K⁺ ions and the carbonyl oxygens (CO) of the filter. The close association of ions K1 and K2 with the carbonyl oxygens of the filter can be seen in Fig. 9. Thus, while the CO(78)-K1 distance increases (i.e., there is loss of a favorable interaction) at ~350 ps, the CO(77)-K1 distance remains almost constant at a value of ~0.27 nm (i.e., the mean value of the K-O distance in a number of small molecule crystal structures; Hille, 1992). Shortly after this transition, the CO(76)-K2 distance increases, as K2 moves from site S3 to S4. Similarly to CO(77)-K1, the CO(75)-K2 distance remains close to its optimal value during this transition. This remains the case until ~930 ps, when K2 leaves site S4 to be replaced by a water molecule. In combination with the structures shown in Fig. 8, these distances suggest that the rings of carbonyl oxygens of residues G77 and T75 distort slightly to track the single file motions of the ions.

View larger version (42K): [in this window] [in a new window]   FIGURE 9   Distances of K+ ions from backbone carbonyl oxygens of the selectivity filter residues as a function of time in simulation MDK3. The vertical arrows indicate the concerted transitions discussed in the text.

It is not only the backbone carbonyl oxygen atoms that play a role in solvating K⁺ ions in the selectivity filter. The O atoms of the T75 side chains at the bottom of the filter, adjacent to the cavity, play an important role in helping to form site S4. As can be seen in Fig. 10, top, when the concerted motion of the K⁺ ions and water takes place at ~350-400 ps, the distance from all four O(T75) atoms to K2 drops to ~0.29 nm, indicating that the ion sits within a ring of threonine side chains. There is also a change in the conformation of the T75 side chains. Their H atoms are initially directed toward the pore to interact with the oxygen atom of a water at or close to site S4 (Fig. 8; t = 100 ps), but they rotate away from the pore when a K⁺ ion occupies the same site (Fig. 8; t = 850 ps). Such a conformational change of channel-lining threonine side chains in response to a permeant cation had been suggested by earlier theoretical studies (Sansom, 1992), and is also seen for the ion in the cavity (see below).

View larger version (35K): [in this window] [in a new window]   FIGURE 10   Distances of threonine side chain hydroxyl oxygen atoms (O) from K+ ions in simulation MDK3. The upper figure shows the development of the interaction of K2 with the ring of T75 side chains that form site S4. The lower figure shows the evolution of the interaction of K3 with two T107 side chains within the cavity. The vertical arrows indicate the start and end of K/O interactions (see text).

Closer examination of the trajectories of water and K⁺ ions along the pore axis for simulation MDK3 (Fig. 11) reveals further changes at the selectivity filter that appear to be coupled to the exit of ion K3 from the channel at ~900 ps (see below). Exit of K3 appears to result in a perturbation of K2 such that it leaves site S4 and moves into the upper end of the cavity. It is replaced by a water molecule (W3 in Fig. 11), which then occupies site S4. Thus, at 1100 ps the configuration within the selectivity filter is water-K1-water-water (for sites S1-S2-S3-S4, respectively). This might suggest that repulsive interactions between K3 and K2 may hold K2 in the filter until K3 leaves the cavity. However, it must be remembered that 1) the interaction between K2 and K3 will be screened by water; and 2) the use of a 1.7-nm cutoff in the electrostatic energy calculations will artefactually remove K2-K3 interactions once K3 leaves the channel. Furthermore, the exit of K3 has only been seen in one realization of the process. Further simulations and a more detailed energetic analysis will be required to resolve this (Shrivastava and Sansom, unpublished results).

View larger version (36K): [in this window] [in a new window]   FIGURE 11   K+ ions and water molecules in the selectivity filter in simulation MDK3, showing locations on the pore axis (z) of K1 and K2 and of three water molecules (of which W2 is the crystallographic water). The horizontal broken lines indicate the approximate extent of the pore, from the intracellular (z = 1.6 nm) to extracellular (z = 5.8 nm) mouth. Arrow a indicates the concerted motion of K1-W1-K2 down the filter; arrow b indicates the exit of K3 from the channel; and arrow c indicates when K2 leaves site S4 of the filter to enter the cavity.

K3 in the cavity

The ion K3 was placed at the geometric center of the pore cavity in the beginning of the simulation (Fig. 1). This is close to the optimal location for K3 discussed by Roux and MacKinnon (1999). However, at ~200 ps K3 moves a little deeper into the cavity. Visualization of the change reveals that the K3 ion moves toward the ring of T107 side chains, and forms a strong interaction with the O atom of two of these (K-O distance ~0.29 nm, compared to a mean distance for K-O in a number of inorganic crystal structures of ~0.27 nm; Hille, 1992). Note that this does not correspond to a large lateral displacement of K3 from the central pore axis. It remains in contact with these two T107 side chains until ~900 ps (Fig. 10, bottom). Thus, threonine side chains also provide a K⁺ ion site within the otherwise largely hydrophobic cavity. At 900 ps (arrow b in Fig. 11), the K3 ion leaves the cavity (and the channel), exiting through the intracellular mouth.

Examination of those water molecules interacting most closely with K3 (Fig. 12) reveals that from ~300 ps onward it forms close interactions with the oxygens of two waters that persist until ~950 ps. Thus, from ~300 to 900 ps the K3 ion interacts with four oxygen atoms, two from T107 side chains and two from water molecules. As it starts to leave the cavity, the interactions with the two T107 side chains are lost, and just before it exits it interacts with four water oxygens (arrow in Fig. 12). All four of these interactions are broken as the ion leaves the channel, to be replaced by solvation with waters outside the channel as it moves into the bulk solvent. Thus, the water molecules solvating K3 once it has left the cavity are not the same ones that solvate it within the cavity. Once K3 has left the cavity it not longer associates with the channel, but appears to remain close to the surface of the bilayer.

View larger version (31K): [in this window] [in a new window]   FIGURE 12   Ion K3 and water molecules in simulation MDK3, showing the distance of K3 from the oxygen atoms of the four water molecules with which it liganded (as indicated by the vertical arrow) just before its exit from the channel cavity via the intracellular "gate" at ~900 ps.

Changes in pore radius associated with ion exit

As mentioned above, at ~900 ps ion K3 leaves its site interacting with two T107 side chains, moves down to the putative gate (corresponding to the V115 side chain ring), and then leaves the channel. By 1000 ps it has entered the aqueous phase at the intracellular face of the bilayer (Fig. 1 B). The exit of K3 from the channel is linked to breathing motions of the protein. In the crystal structure the intracellular mouth of the pore is narrower than the radius of a K⁺ ion. The pore radius profile, calculated at 5-ps intervals from 475 to 500 ps for MDK3 (Fig. 13 A), illustrates this occlusion of the pore, which is maintained throughout most of this simulation. In particular, the V115 side chain ring maintains a narrow and hydrophobic mouth to the channel. Indeed, short (~5-ps) time scale fluctuations in this region can reduce the pore radius to <0.1 nm. Similar fluctuations in the pore radius profile in the gate region are seen in simulation MDK0 (Fig. 13 B), suggesting that this property of the channel is independent of the presence/absence of K⁺ ions. The pore radius profile of MDK3 also reveals two narrow regions in the selectivity filter. This suggests that small dynamic changes in conformation must occur as the ions and water migrate between sites in the filter, as can be seen from careful comparison of the 100 ps vs. 850 ps structures in Fig. 8 or by visualization of the dynamics of the filter region (Shrivastava et al., 1999b).

View larger version (38K): [in this window] [in a new window]   FIGURE 13   Pore radius profiles sampled at times 475-500 ps (A, B) and 925-950 ps (C, D) for simulations MDK3 (A, C) and MDK0 (B, D). Each graph contains six lines, corresponding to profiles calculated from structures saved at 5-ps intervals. The extent of the channel molecule is indicated by the horizontal gray bar running from z = 1.6 nm to z = 5.8 nm. The position of the V115 side-chain ring (which constitutes the "gate" at the intracellular mouth) is indicated by the vertical arrow. The dotted horizontal lines indicate the Pauling radius of a K+ ion (i.e., 0.13 nm).

In simulation MDK3 the radius profile around the gate changes at the time of exit of the ion (~925-950 ps; Fig. 12 C). The gate opens, as revealed by a 0.03 nm increase in radius to ~0.16 nm, i.e., greater than that of a K⁺ ion. This opening is only seen while K3 leaves the pore, and is maintained for ~100 ps. Analysis of fluctuations in the protein structure suggests that this opening is related to a small movement of the M2 helix of subunit B of the channel. The C RMS fluctuations of this helix, averaged over the entire 1 ns, are ~0.07 nm, compared to ~0.04 nm for the other three M2 helices. Thus small fluctuations in packing of the M2 helices during the course of the simulation may transiently open the gate at the V115 ring that otherwise occludes the pore. In contrast, in MDK0 (Fig. 12 D) the radius remains <0.13 nm at both the intracellular and the extracellular mouths throughout the simulation.

Detailed examination of the configuration around K3 as it is leaving the channel shows that it sits in a hydrophobic pocket. It interacts above with two waters in the cavity (both K-O distances = 0.27 nm) and below with two waters in the intracellular mouth (K-O distances = 0.28 and 0.29 nm). However, it does not appear to form any favorable interactions with polar groups of the protein, either main chain or side chain. Thus, even while the gate is open there is likely to be an energetic barrier to ion translocation. Inspection of the potential energies (not shown) of K3/protein and K3/water interactions reveal a loss of favorable electrostatic interactions between ion and protein as K3 passes through the gate, which is more than compensated for by a gain in favorable interactions with water. However, a proper analysis of energetics would require estimation of free energies (Roux, 1996; Roux and Karplus, 1991).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS CONCLUSIONS REFERENCES

Biological relevance

Overall, our simulations reveal that interactions of K⁺ ions and water with the KcsA channel at both the selectivity filter and at the intracellular gate are dynamic. At the selectivity filter a concerted translocation of K1-W2-K2 between sites takes place. This is associated with small "peristaltic" changes in the conformation of the filter backbone to maintain optimal K-O interactions (as was previously seen in simulations of gramicidin-cation interactions; Roux and Karplus, 1991). Ion-ion interactions within the filter have been analyzed by Dorman et al. (1999), who also commented on the inherent flexibility of the structure. The MDK3 simulation confirms that there is a destabilizing interaction between adjacent K⁺ ions within the pore (data not shown) as originally suggested by Doyle et al. It appears that long-range interactions may occur between the ions in the filter and that in the cavity, although interpretation of this is complicated by use of a 1.7-nm cutoff in the electrostatic energy evaluation. To some extent these destabilizing interactions may be overcome by the favorable interaction of K3 with two T107 side chains. However, the movement of K2 toward the cavity after K3 has left the channel suggests that movements in the selectivity filter and cavity may be coupled. Although a free energy analysis is needed to probe these effects in more detail, our simulations clearly reveal the dynamic and concerted behavior of K⁺ ions and water within a multi-ion pore. Furthermore, recent simulations on a homology model of an inward rectifier potassium channel (Capener, Shrivastava, and Sansom; unpublished data) reveal similar concerted K⁺/water movements to those seen in the current study. This suggests that this aspect of our simulation results is robust to changes in the detail of the channel structure, and so may be a more general property of potassium channels.

The concerted, single-file motion within the selectivity filter is consistent with both structural and electrophysiological data. In particular, the x-ray studies (Doyle et al., 1998) data indicate partial occupancy of sites S3 and S4 by a cation, demonstrating that multiple patterns of ion/filter interactions are possible, even within a crystal. Potassium channels have long been suggested to form single-file, multi-ion pores (reviewed by Hille, 1992). This is consistent with the observations of ion movement through the selectivity filter in the current simulation. Furthermore, although in principle the cavity might be occupied by multiple ions (but cf. the electrostatics calculations of Roux and Mac-Kinnon, (1999) which would argue against this), the mutual repulsive interactions of K2 and K3 may suggest that the cavity essentially acts as an extension of the single file pore. We note (as have several others, e.g., Wallace, 1999) the resemblance between ion/water movement in the KcsA selectivity filter, and the concerted single-file motions of, e.g., water (Chiu et al., 1999) in gramicidin channels.

The events at the gate in simulation MDK3 are of interest in the context of spin-labeling studies (Perozo et al., 1998, 1999) which suggest that activation of KcsA at low pH is linked to a widening of the pore in this region. Our results suggest that relatively small (~0.1 nm or less) changes in the region of the V115 ring can open the channel. The simulation results correlate nicely with the spin label data (Perozo et al., 1999) in that the latter reveal a maximal change in the vicinity of residues 116 and 117, i.e., next to the V115 "gate." Of course, it is possible that larger (and presumably slower) changes may occur during physiological gating. Such changes have been suggested by Perozo et al. (1999). Indeed, for Kv channels the results of Armstrong (1971) and Liu et al. (1997) suggest that the gate may open more widely than in the present simulation, providing relatively unhindered access to the cavity for quaternary ammonium ions or cysteine-directed reagents, respectively. Interestingly, the S6 helix of Kv channels (which corresponds to the M2 helix of KcsA) contains a highly conserved PVP sequence motif. Simulation studies on single S6 helices (both in vacuo (Kerr et al., 1996) and inserted into a lipid bilayer (Shrivastava et al., 1999a)) suggest that the PVP motif may provide a molecular hinge. This might provide a mechanism for opening of the gate of Kv channels wider than is seen in the current simulations.

On the basis of simulations of KcsA embedded in a slab of octane, Guidoni et al. (1999) noted the importance of a salt bridge between D80 and R89 of neighboring subunits in stabilizing the tetrameric structure of KcsA. Although not present in the x-ray structure, we note that this salt bridge formed, at all four subunit interfaces, in MDK3 during the early stages of the simulation and was maintained throughout. This merits further investigation, in particular with respect to the role of side chain ionization states on the behavior of the channel.

A further aspect of these simulations merits comparison with both experimental data and simulation studies. In the absence of K⁺ ions, e.g., in simulation MDK0, the selectivity filter undergoes significant conformational changes. It is encouraging that Guidoni et al. (1999) also see enhanced flexibility in the structure of the selectivity filter in the absence of ions, despite a rather different simulation system. It has been suggested that C-type inactivation of potassium channels may be impeded by ions bound near the external mouth of the channel (Ogielska and Aldrich, 1999), and that in the absence of K⁺ ions some potassium channels enter a "defunct" state (Khodakhah et al., 1998; Melishchuk et al., 1998). Our simulations suggest that a possible early event in this process is a change in conformation of the selectivity filter. Of course, subsequent slower conformational transitions (which cannot be captured in our current simulations) may occur which link the changes at the selectivity filter to C-type inactivation per se.

Limitations of methodology

One of the main limitations of our simulations is the use of a cutoff (albeit 1.7 nm) for long-range electrostatic interactions. As discussed above, this cutoff means that in MDK3, the K1 ion does not "see" the K3 ion. It is unclear what the optimal solution is. Either one could increase the cutoff distance to, say, 3.0 n, or one might use a more sophisticated approach to long-range electrostatics, such as Ewald summation (Tieleman et al., 1997; Tobias et al., 1997). However, as the lipid parameters used have not been evaluated using Ewald summation, a series of detailed control simulations will be required to identify an optimal protocol. Until such simulations are available, one must retain a degree of caution in interpretation of the current results. A second limitation is that the simulation system (area ~9.5 nm × 9.5 nm) is still relatively small. As emphasized by Gouliaev and Nagle (1998), on a larger (50-100 nm) scale, distortions (undulations, compression/expansion) of the bilayer may occur. This suggests that channel gating might be stochastically modulated as a protein travels between different regions. At present it would be too challenging to simulate a sufficiently large system to allow this to be explored directly by MD simulation. However, simulations with different lipid species could reveal the extent to which the channel protein dynamics are sensitive to the local bilayer environment.

One limitation that may be addressed by further simulations is the relatively short time scale (1 ns). Future work will extend the current simulations by an order of magnitude (Shrivastava and Sansom, work in progress). This may reveal further aspects of M2 helix movement in relationship to channel gating. In particular, longer simulations may provide some clues as to the rigid body motions of the M2 helices and their relationship to gating. Longer simulations will also enable detailed exploration of whether the protocol adopted here has allowed the bilayer to fully re-equilibrate after "insertion" of the protein.

Future directions

This study has not addressed the physical basis of the K⁺ ion selectivity of KcsA. However, given the deformability of the selectivity filter in our simulations it seems likely that a straightforward stereochemical explanation of K⁺ vs. Na⁺ selectivity may prove too simple. An alternative approach may be to calculate free energy profiles for K⁺ vs. Na⁺ ions moving through the KcsA channel. However, this will present a number of methodological challenges, particularly for a multi-ion channel.

So far we have omitted any representation of a transbilayer voltage from our simulations. It will be of interest to repeat our simulations with a difference in voltage across the membrane to see whether ion movement through the channel may occur in either direction (as would be expected from a physiological standpoint). Such simulations have been performed for bacterial porins (Suenaga et al., 1998). Simulations on simpler channel systems (Zhong et al., 1998c) have suggested how a transbilayer voltage might be included in bilayer simulations. It will also be important to increase the number of ions in the bulk aqueous phases on either side of the bilayer. Furthermore, one must remember that 22 residues are absent from the N-terminus and 35 from the C-terminus of the KcsA chains in the x-ray structure. It would be interesting to attempt to include these in future simulations, although Perozo et al. (1999) have shown that C-terminal truncation of KcsA does not lead to loss of gating. Despite these reservations, it is encouraging that the current simulations have provided direct visualization of fundamental events of underlying ion permeation through a potassium channel in a phospholipid membrane.