INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS REFERENCES

The bacterial potassium channel KcsA marks the first permeation pathway high-resolution structure (Doyle et al., 1998) from the superfamily of the voltage-gated ion channels, which may be present in most living cells and underlie innumerable excitability, transport, signaling, and osmoregulatory functions (Hille, 1992). The primary open-closed gating of KcsA is by pH lowering (Heginbotham et al., 1999) with a voltage-dependent distribution of conducting open substates (Meuser et al., 1999). Perozo et al. (1999) used electron paramagnetic resonance spectroscopy (EPR) with site-directed spin labeling to show that during gating the transmembrane helix, which forms the lumenal surface of the pore, undergoes a rotational motion that provides a wider opening at the intracellular end of the channel.

A number of computational studies using methods of molecular dynamics (MD) (Allen et al., 1999, 2000; Åqvist and Luzhkov, 2000; Bernèche and Roux, 2000; Shrivastava et al., 2000; Shrivastava and Sansom, 2000), Brownian dynamics (BD) (Chung et al., 1999; Corry et al., 2000), and electrostatics methodologies (Roux and MacKinnon, 1999; Guidoni et al., 2000) have been done on KcsA, providing valuable information about the biophysics of the channel. General features of permeation as shown by these studies include the stabilization of the selectivity filter by ion occupancy, the obligatory single filing of ions through the selectivity filter, and the preferential attraction of cations into the selectivity filter by a combination of side chain and backbone charges exposed to the lumen. It has also been calculated that ions are significantly attracted into the wide cavity in the center of the channel by the dipole of the pore helix (Roux and MacKinnon, 1999).

However there has not yet been a quantitatively successful description of ion permeation in this channel that predicts electrophysiological properties by applying physics with no arbitrary parameters directly to experimentally measured structural data. The purpose of the present study is to describe the implementation of a hierarchical strategy (Jakobbson, 1998) combining multiple computational methods in a coordinated way in which all the parameters in the calculations are derived directly from structural data with no arbitrarily adjustable parameters and to apply this strategy to the KcsA potassium channel. As a high-resolution open KcsA channel structure is not yet available, we have constructed a series of structures of different opened sizes based on our best interpretation of published experimental data for how the open channel varies from the x-ray structural model, which is functionally closed. The resulting computed ion fluxes accurately reproduce known electrophysiologically measurable properties for this channel.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS REFERENCES

Molecular modeling of the open channel

Our strategy for molecular modeling of the channel takes several factors into account. 1) A motion of the M1 transmembrane helix and pore-lining M2 transmembrane helix tending to open the intracellular end of the channel, as suggested by the EPR spectroscopy of Perozo et al. (1999). 2) A degree of opening of the intracellular end of the channel sufficient to permit tetrabutylammonium ion (TBA) to be an optimal channel blocker relative to other quaternary amines (Meuser et al., 1999; Splitt et al., 2000) and to enter the channel completely into the deep cavity (Zhou et al., 2001). 3) The preservation of the shape of the extracellular vestibule, as suggested by experimental data, that differences in toxin binding affinity of open and closed K channels is a function of interactions among K ions in the channel rather than a change in the shape of the vestibule (Terlau et al., 1999). 4) The maintenance of hydrophobic matching of the transmembrane helices with the surrounding implicit membrane, by maintaining a constant helix tilt with respect to the bilayer normal, which is the same as the channel axis. The tendency to preserve hydrophobic matching is a generally accepted principle of protein-lipid interactions, albeit not always perfectly obeyed (Killian, 1998). The specific modeling process is described below.

The x-ray structure of KcsA in the Protein Data Bank is an L90C mutant with missing residues 1 through 22 and 120 through 158 and truncated side chains at Arg-27, Ile-60, Arg-64, Glu-71, and Arg-117 (Doyle et al., 1998). This structure was refined by completing the truncated side chains and optimizing their positions using the program Modeler (Sali et al., 1990). The ends of the four subunits were capped with the nontitratable groups CH₃ ---C==O and NH---CH₃ as extensions of the -helices. A series of energy minimizations was done using a distance-dependent dielectric, leading to a structure that satisfied the structural criteria of Procheck (Laskowski et al., 1993) and information-based probability density functions (Wall et al., 1999).

To achieve an opening of the intracellular end of the channel consistent with the hydrophobic matching and maintenance of the external vestibule shape, the extracellular ends of the transmembrane helices M1 and M2 were kept fixed, and the M2 helices were rotated about an axis parallel to the membrane normal (channel axis). The M1 helices were rotated in the same fashion as necessary to eliminate steric overlap with M2. The M1 was rotated approximately the same amount as was the M2.

The refined crystal structure and a structure formed by a 20° rotation of the M2 were inspected for gaps in packing as visualized by RasMol (Sayle and Milner-White, 1995) in space-filling mode, as seen in Fig. 1, A-D. When viewed from the extracellular end and from within the membrane plane, both the closed and opened structures have some gaps between the M1 and M2, which are presumably filled by hydrocarbon when the channel is in a membrane. Our construction procedure may widen existing gaps but introduces no substantial new gaps into the protein structure.

View larger version (85K): [in this window] [in a new window]   FIGURE 1   Packing of transmembrane helices in the closed (A and C) and modeled 20°-opened (B and D) KcsA structures (M2 helices, dark gray; M1 helices, light gray). The P regions of the channel are omitted for clarity. Views in A and B are from the extracellular end of the channel. Views in C and D are side-on from the membrane plane such that the interhelical gaps appeared to be maximized. (E) Channel pore diameter given by the protein van der Waals surface versus position along the channel axis, as measured by the program HOLE (Smart et al., 1993), for several M2 helix rotations (x-ray model, 0°). The protein extends from 0-60 Å with the extracellular end oriented at left. The narrow region from ~10-25 Å contains the selectivity filter, and the narrow region from ~30-50 Å in the x-ray structure is lined with aliphatic side chains from the four M2 helices. (F) Superimposed ribbon representations of the channel in the closed (gray) and 20°-open conformation (dark gray) as viewed from the membrane plane (spheres, light gray) for electrostatics calculations. The remaining sections corresponding to the selectivity filter, extended, and turret regions (black) are common to both structures. Orientations in panels C, D, and F are with the extracellular end at top.

Fig. 1 E shows the diameter of the permeation pathway as determined by the program HOLE (Smart et al., 1993) and its variation along the axis in the Doyle et al. (1998) x-ray structure and in structures of various opened sizes. (For convenience, the term "opened channel" refers to the 20°-opened structure unless otherwise stated. The 20° estimate for the degree of rotation of the M2 helix in the fully opened channel is based on the estimate from Perozo et al., 1999.) Two very narrow regions in the x-ray structure are found at the selectivity filter and near the intracellular end of the channel. The intracellular end widens to an average diameter of ~13.5 Å in the opened channel. This size is consistent with experimental evidence on channel blocking by quaternary ammonium ions (Meuser et al., 1999; Splitt et al., 2000), which suggests the size of the intracellular opening for the maximally conducting substate of the channel is larger than TBA. Further evidence that the channel opens widely enough to permit TBA to enter comes from a recent x-ray structure of the channel-TBA complex showing that TBA enters completely into the channel (Zhou et al., 2001). (The K channel structure complexed with TBA is not much different from the closed structure, but a transient opening at the intracellular end would be required for the TBA to enter.) Estimates for the effective molecular diameter of TBA vary from a minimum of ~9.5 Å (Splitt et al., 2000) to a value of 11.6 Å based on an idealized geometry (Huang et al., 2000).

It should be noted that our open channel structure is not certain in complete detail. However, there is strong experimental evidence that the intracellular end of the channel opens to approximately the extent that we suggest, whereas the extracellular vestibule and selectivity filter do not change in dimension or conformation during channel opening. On the other hand, the M2 residues are generally farther away from the channel axis in our opened channels as compared with the closed channel, whereas EPR experiments (Perozo et al., 1999) suggest that some of the residues move closer to the axis. This might be achieved by adding to the motion that we have suggested a rotation of the helix about its axis and/or a kink within the helix itself. Given the absence of direct evidence for specific movements of these kinds and the subtleties inherent in interpretation of EPR data, we judged it not useful to pursue further structural modeling without more detailed structural data on the open state. Because of the residual uncertainty of the structure, as well as other approximations in the calculations that we discuss later in the paper, we emphasize that our results should be considered as only semiquantitative representations of the channel behavior. We believe the major conclusions of the paper depend only on the intracellular end of the channel opening to approximately the same degree as in our model with the selectivity filter unchanged and will therefore stand as long as the correct model of the open channel has these features.

Electrostatics: pK_a calculations and potential profile for ion permeation

Electrostatics calculations were used to calculate the ionization states of protein residues under various physiological conditions. We used the method described in Gilson (1993) and Antosiewicz et al. (1994) as implemented in the University of Houston Brownian Dynamics (UHBD) suite of programs (Madura et al., 1995) used to calculate electrostatic energies. UHBD requires the specification of dielectric constants for the protein and surrounding electrolyte. To emulate the low-dielectric environment around the channels for these calculations, a layer (thickness, 23.4 Å) of neutral nonpolar spheres (radius, 2.0 Å), postulated to be situated in between the aromatic residues at the lipid-water interfaces (Cowan et al., 1992), was placed around each structure, as shown in Fig. 1 F. The spheres were packed at hydrocarbon densities out to ~37 Å from the channel axis. The protein dielectric constant was set to 20, which appears to be the optimum for computing the ionization states of residues in many proteins (Antosiewicz et al., 1994), and the nonpolar spheres were assigned a value of 20 as well. Outside the regions delineated by protein and membrane the dielectric constant was set to 80 to represent water. UHBD automatically accounts for the image charges arising from introduction of a dielectric interface.

The partial charge parameters used in the calculation of electrostatic energies were taken from the CHARMm (Brooks et al., 1983) parameter set and radii from the OPLS (optimized parameters for liquid systems) parameter set (Jorgensen and Tirado-Rives, 1988), following Antosiewicz et al. (1994). Reference pK_a values corresponding to the pK_a values of model compounds in solution for each titrating residue type are taken from Antosiewicz et al. (1994). A modified approach (Tanford and Roxby, 1972) was used to determine effective pK_a values in which an ionization polynomial is evaluated exactly for clusters of residues with significant charge-charge correlations and a mean field approximation is used for less significant intercluster interactions (Gibas et al., 1997). The cluster method requires that each histidine be assigned a tautomer, and in all cases the tautomer in which N2 is the titrating site was used. The calculations used a finite-difference method for computing electrostatic potentials at the lattice points of a grid of spacing 1.2 Å surrounding the titratable site. Focusing grids of sizes 1.0, 0.75, and 0.25 Å were used to refine the potentials. (The position and orientation of the lattices were fixed with respect to the coordinate axes for all channels studied.) The pK_a values were computed using the linear Poisson-Boltzmann electrostatics module within UHBD. Test calculations with the nonlinear Poisson-Boltzmann method revealed that the modifications of the ionization states were not large enough to justify their greater computational demand. The ionization states were not substantially affected by either the presence of a single on-axis K⁺ ion or by ionic strength in the range 0.15 to 1.50 M. Symmetric pH conditions were used, and all calculations were carried out at 298 K.

For computing potential profiles of a single ion along the channel axis, the calculated ionization states were converted into new, effective partial charges for the side chain atoms. The total force acting on the ion as a function of position along the axis was then calculated using UHBD. The axial component of this force was integrated over the length of the channel to give the potential profile of the ion. Because of the asymmetry of charge in the channel protein combined with the absence of shielding charges in the electrolyte in this reduced system, it was necessary to normalize the computed transmembrane potential to zero by the addition of a constant field.

Molecular dynamics: determining the mobility of ions within the channel

MD simulations of ion motions throughout the channel were performed using the program GROMOS96 (van Gunsteren et al., 1996) with solvation of the refined version of the channel by methods we previously used for the gramicidin channel (Chiu et al., 1993) and for part of the permeation pathway for a model of the sodium channel (Singh et al., 1996). Specifically, the wide parts of the channel pore were hydrated, and the ends of the channel were bathed in approximately hemispherical caps of water extending to the aromatic collars. Three potassium ions, two of which were initially located in the outer crystallographic positions at either end of the selectivity filter and one ion located in the central cavity, were included for equilibration for 1 ns at 298 K. To prepare the system for ion diffusion measurements, two of the ions were replaced with water molecules, and the system was re-equilibrated while strongly restraining the restraining ion to the channel axis. The position of the strong restraints was updated in 1-Å increments along the channel axis, allowing for reequilibration in each instance. Several reequilibrated systems were selected to represent a subset of ion positions along the channel axis, including those in which the ion was located in the outer vestibule, within the selectivity filter, in the central cavity, and in the narrow, hydrophobic part of the channel. Production runs of 27 ps with no restraints on the ion were generated in each region of the channel.

Velocity autocorrelation functions were calculated for the unrestrained ion trajectories in each part of the channel. The velocity autocorrelation function (McQuarrie, 1976) relevant to one-dimensional motion of ions along the channel axis (see below) is given by the expression v_z(t₁)v_z(t₂), in which v_z is the z-component of the velocity of the ion at two times, t₁ and t₂, and the angle brackets denote an equilibrium average. In an equilibrium system the correlations depend on only differences in time t = t₂  t₁. We fit the velocity autocorrelation function to the single exponential function

(1)

in which A₀ is a scaling factor and the time constant, , is related to the fluid dynamic diffusion coefficient, D_f, by the expression

(2)

in which m is the mass of the ion, T is the temperature, and k_B is Boltzmann's constant.

Brownian dynamics: calculating fluxes and distributions within the channel

BD simulations were performed using methods we have previously reported (Cooper et al., 1985; Bek and Jakobsson, 1994) to compute ion fluxes. These calculations are one-dimensional, in effect constraining the ions to move along the channel axis. For each time step, each ion in the system is moved a certain distance based on a combination of deterministic motion due to a free energy gradient and random motion representing thermal fluctuations, according to the equation

(3)

in which z is the distance that an ion moves in the time interval t, D is the ion diffusion coefficient, R is the ideal gas constant, T is the absolute temperature, U/z is the free energy gradient, and  is a random number chosen from a normal Gaussian distribution centered at zero. The free energy gradient, U/z, is approximated by its electrostatic component z_ionF(V/z), in which z_ion is the valence of the ion, F is Faraday's constant, and V/z is the voltage gradient. The voltage gradient is the force seen by an ion due to the channel and surrounding solvent as it moves through the channel, plus the applied transmembrane voltage, plus the electric field due to other ions in the channel.

A summary of the parameters for the BD calculations is listed in Table 1. Based on the MD results (described in the Results section below), we set the K⁺ diffusion coefficient inside and outside the channel equal to the dynamic ion diffusion coefficient for K⁺ in bulk water. A dielectric constant of 20 was used for ion-ion interactions, which is consistent with the ions being in a narrow cavity surrounded by protein. The time step was 0.4 ps, and the time step and mobility together lead to a mean thermal jump distance of 0.45 Å per time step. Ions enter the channel by means of the "entrance tube" algorithm (Cooper et al., 1985; Jakobsson and Chiu, 1987; Bek and Jakobsson, 1994), a particular implementation of the flux-over-population method (Farkas, 1927, as cited by Hänggi et al., 1990). The capture radius for ions attempting to enter the channel was taken as the average radius of the intracellular half of the channel (Fig. 1E). Thus the BD simulations contain no arbitrarily chosen parameters.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS REFERENCES

Computing potential profiles

An overview of our open channel models is shown in Fig. 1, A-F. Detailed information about the generation of these structures is given in the Materials and Methods section. Given these models, we first compute the ionization states of the titratable residues in both the x-ray and opened channel structures at pH 4 and 7, as it has been demonstrated experimentally that the channel is functionally closed at pH 7 and is fully open at pH 4 (Cuello et al., 1998; Heginbotham et al., 1999). The results in Table 2 show that the glutamates, located in the vestibule and toward the C terminus, undergo partial protonation at low pH, whereas the histidines, toward the N terminus, carry a full charge. (For comparison, the values used by the GROMOS96 simulation package standard force field are also given and are in close correspondence with the computed ionization states at neutral pH for all residues except for Glu-71.) Ionization states of side chains at various ionic strengths were computed similarly. We found that neither the channel conformation nor variation of the ionic strength had an appreciable effect on computed ionization states.

Fig. 2 shows the potential profiles as seen by a K⁺ ion traversing the channel along the channels axis. The potential profile for the x-ray conformation at pH 7 has three notable features: a deep potential well containing the selectivity filter, a second potential well near the intracellular end of the channel, and a pronounced potential barrier in the narrow region of the channel. All curves show a deep potential well near the extracellular end of the channel that is due to the negatively charged carbonyl oxygens in the selectivity filter and, to a lesser extent, to the negatively charged residues in the vestibule of the channel. This well is somewhat shallower at pH 4 than at pH 7 due to the partial neutralization of Glu-51 in the vestibule. The potential well near the intracellular end of the x-ray structure is due to negatively charged Glu-118 and shows a large pH dependence. The pH dependence here for the open structure is much less because Glu-118 is now farther from the channel axis. The barrier, resulting from the dielectric contrast between the high-dielectric electrolyte and the low-dielectric protein and surrounding membrane, is an image potential barrier (Sancho et al., 1995) that can be thought of as a Born energy, or desolvation energy, associated with moving an ion from the surrounding electrolyte into the narrow nonpolar channel interior. The profile for the opened channel at pH 4 reveals that this image potential barrier is greatly reduced upon widening the channel lumen. It is seen that the depth of the potential well at the selectivity filter is not very sensitive to the degree to which the channel is opened but is sensitive to the ionization state of the protein. The potential profiles as computed for other ionic strengths (not shown) provide similar results. Fig. 2 B shows the potential profile for the potassium ion traversing the channel axis for different sizes of channel openings at pH 4. It is seen that the potential barrier near the intracellular end of the channel is systematically reduced as the channel is opened further.

View larger version (27K): [in this window] [in a new window]   FIGURE 2   (A) Potential profile of a potassium ion in translocating the channel along the channel's axis in the x-ray and open channel structures at neutral and acidic pH using partial charge data from Table 2. Those profiles in bold relief correspond to physiologically relevant conditions and are the basis for the subsequent calculations. (B) Potential profiles for K+ in the structures of intermediate open sizes (as measured by the degree of M2 rotation) at pH 4. Note the decline in the Born energy barrier near the intracellular end of the channel as the channel is opened.

In summary, our modeling and electrostatics calculations thus show that the depth of the potential well at the selectivity filter is controlled by the ionization states of titratable residues of the protein, whereas the image potential barrier near the intracellular half of the channel is controlled by the size of the pore at the intracellular end.

Computing ion mobilities

Velocity autocorrelation functions were calculated from the ion trajectories generated by MD simulations. Values for the fluid dynamic diffusion coefficient (due to local friction), determined from Eqs. 1 and 2, were found to range from 1.8 × 10⁵ to 3.3 × 10⁵ cm² s¹, in close agreement with the experimental value of 2.5 × 10⁵ cm² s¹ that pertains to K⁺ ions in bulk water (Hille, 1992). Because of the complexity of the ion solvation environment, the velocity autocorrelation functions were more complex than simple exponentials. However, the best fit to an exponential still suffices to give an approximately correct diffusion coefficient. In the wide part of the channel, the ion is in fact solvated by water. In the narrow part of the channel, it appears from our results that the carbonyl oxygens substitute approximately equivalently to bulk water oxygens in determining the local friction for ion movement. It should be pointed out that the diffusion coefficient in this case is not the effective diffusion coefficient, but rather the reciprocal of the local friction. The effective diffusion coefficient in the selectivity filter is certainly less than in bulk water, due to the potential well created by the carbonyl oxygens, even while the local friction is the same as bulk water. This situation is similar to that pertaining in gramicidin, where we earlier found that the local friction for sodium ions in gramicidin was similar to bulk water, whereas the effective diffusion coefficient was lower by a factor of 10, due to the potential wells created by carbonyl oxygens in that structure (Chiu et al., 1993).

Computing ionic fluxes

Ionic fluxes were calculated using a one-dimensional BD method for the ions. The resulting conductance calculated for the closed channel at both pH 7 and 4 was a fraction of a picosiemen. In Fig. 3 A the computed current-voltage (I-V) curve for the fully opened channel for 250 mM K⁺ shows an approximately linear relationship over the voltage range of ±100 mV in agreement with the experiment. The corresponding average channel conductance over this voltage range is ~110 pS, in reasonable agreement with experimentally determined maximal conductances of 90 pS (Schrempf et al., 1995) and 135 pS at 200 mM K⁺ (Cuello et al., 1998) and of ~130 pS at 250 mM K⁺ (Meuser et al., 1999). The mean number of ions in the channel and the Ussing flux ratio exponent were computed for various transmembrane voltages (Table 3) and were found to be close in value to each other, in agreement with the original theory relating these two quantities (Hodgkin and Keynes, 1955). Although there are currently no experimental flux ratio exponent data for KcsA, a recent study (Stampe et al., 1998) on an inward-rectifying K channel has shown ratios between 2.1 and 2.5 for membrane potentials ranging from 50 to 25 mV, which is in close agreement with our simulations.

View larger version (18K): [in this window] [in a new window]   FIGURE 3   (A) Current-voltage curves as computed by BD for the opened channel at pH 4 for 250 mM K+ using the parameters in Table 1, compared with experimental results (Meuser et al., 1999) for the largest conducting substate under comparable conditions. Each simulated point represents a 0.2-ms-long simulation. (B) Current-concentration curves for the opened state at pH 4 for several membrane potentials with parameters from Table 1, compared with experiment (Meuser et al., 1999). Each curve is normalized to the respective value of the current at 250 mM K+. (C) Current-voltage curves as computed by BD for channels of different opened sizes at 250 mM K+. Two factors are responsible for the difference in currents. One is the size of the Born potential barrier near the intracellular end of the channel, as shown in Fig. 2 B. The other is the capture diameter for ions to impinge on the intracellular end of the channel, as visualized in Fig. 1 E. The capture diameters used for the calculations shown were 5.5 Å for 5°, 8.4 Å for 10°, 11.0 Å for 15°, and 13.5 Å for 20°. In all calculations shown in Fig. 3 the standard deviations of the simulated currents are smaller than the size of the data symbols used.

Fig. 3 B shows the relationship between current (normalized to the 250 mM K⁺ value) and concentration for the opened channel at several transmembrane voltages with comparison with experimental values (Meuser et al., 1999). (The channel-ion potential profiles were recalculated for each ion concentration.) Whereas the calculated concentration dependence of the simulated current is similar to that seen in experiment, the simulated currents tend to saturate at a slightly lower concentrations. We believe this behavior may be due to the approximation of one-dimensional ion motion implicitly used in the BD simulations, an approximation that is not as good at near-saturating concentrations as it is at lower more physiological concentrations.

Fig. 3 C shows the I-V curves for single channels for the different open sizes shown in Fig. 1 E. The modeled partially opened channels produce intermediate conductances in a fashion similar to conducting substates. The conductance increase as the pore size increases is partly due to the reduction in the Born potential barrier near the intracellular end (Fig. 2 B) and partly due to the increase in the capture diameter of the channel when the channel opening is larger. It may be that the basis of conducting substates in K channels is the size of the pore aperture at the intracellular end of the channel in a similar fashion to this model.

Features of Brownian ion trajectories

The sample of BD ion trajectories in Fig. 4 demonstrates the extreme ability of the negative charges lining the selectivity filter to concentrate cations. Here the ribbon structure of the opened channel is superposed on the trajectories so that the structural basis for the preferred ionic positions can be seen. There is a clear preference for a double-occupancy state in the selectivity filter, consistent with the computed mean number of ions in the channel (Table 3). The histogram of ion positions computed for the entire length of the BD simulation reveals two prominent peaks in the selectivity filter separated by ~7 Å, similar to the pattern of ion distribution as revealed by x-ray crystallography (Doyle et al., 1998). Fig. 4 also shows the occasional appearance of a single-occupancy state in the selectivity filter where the mean position of the single ion lies in the carbonyl "cage" in between the two that are preferred during double occupancy. This relationship between singly and doubly occupied states is similar to the findings of recent free-energy perturbation MD simulations (Åqvist and Luzhkov, 2000).

View larger version (45K): [in this window] [in a new window]   FIGURE 4   A 120-ns sample of BD ion trajectories in the opened channel at pH 4 at 0 mV membrane potential in 250 mM K+ solution using parameters in Table 1. For ease of visualization the data were smoothed using a window of 120 ps, effectively filtering out oscillations greater than ~1010 Hz. Ions that have successfully traversed the channel in this time window are shown in bold relief. A RasMol (Sayle and Milner-White, 1995) ribbon representation of the backbones (gray) from two opposing subunits showing the negatively charged carbonyl groups (black) projecting into the selectivity filter lumen is superposed.

The trajectories in Fig. 4 also reveal that two ions occupying the selectivity filter move nearly in unison in response to an incoming third ion. The movement of one ion out of one end of the selectivity filter and the entrance of another from the other end is so close to simultaneous that on the time scale of Fig. 4, the trajectories suggest a "knock-off" mechanism. Recent experimental data on channel blocking of KcsA by barium ions has suggested similarly that a third ion could enter the selectivity filter to push the queue of ions along (Jiang and MacKinnon, 2000).

Because the selectivity filter holds two ions practically all the time, and because the mean channel occupancy is only fractionally over two (Table 3), it follows that the wide part of the channel hardly ever holds more than one ion. This shows why the one-dimensional approximation to the motion is reasonably accurate. In the selectivity filter the motion really is one-dimensional because the channel is narrow there. In the wide part of the channel, only the projection of the motion of the single ion onto the channel axis is all that contributes to the flux. Thus, ignoring the lateral motion of the ions in this region does not introduce gross inaccuracies into the calculation.

Ion channels are essentially electrical resistance elements (Hodgkin and Huxley, 1952). In the KcsA (and presumably other) potassium channels, the selectivity filter is in series with the wider portions of the channel at the extracellular and intracellular ends. We analyzed this channel as two resistors in series, one comprising the selectivity filter and the other the rest of the channel. Because the selectivity filter is almost always doubly occupied, regardless of bath electrolyte concentration, we postulate that the resistance of the selectivity filter is a constant. We further postulate that the resistance of the rest of the channel is inversely proportional to electrolyte concentration, similar to bulk electrolyte. This leads to the functional form

(4)

in which R is the resistance of the single channel (reciprocal of the conductance), R_f is the intrinsic resistance of the selectivity filter, and A is a constant. This model of resistors in series provides a physical basis for the frequently observed property of ion channels to follow Michaelis-Menten kinetics (Aidley and Stanfield, 1996). The solid and dashed lines in Fig. 5 show the fit of Eq. 4 to both computed and experimental data (same data as Fig. 3 B) for a voltage of 100 mV. It is seen that the fit is very good. Similarly good fits to Eq. 4 are obtained for the full range of simulated voltages. In this view, the reciprocal of R_f (~500 pS) is the theoretical maximal conductance achievable by a potassium channel if all resistance in series with the selectivity filter could be eliminated. The ratio of selectivity filter resistance to that of the wider pore regions was computed (Fig. 5, inset). At low bath potassium concentrations, the majority of the KcsA resistance appears to be in the wider regions of the channel, whereas at high concentrations the majority of the resistance is in the selectivity filter.

View larger version (25K): [in this window] [in a new window]   FIGURE 5   Opened-channel resistance, R, computed from the ratio of voltage to current, as a function of electrolyte concentration at 100 mV for BD simulations (circles) and for electrophysiological measurements of the maximal conducting substate (Fig. 4 of Meuser et al., 1999) (triangles). Fitting of the data sets to Eq. 4 yielded constants of Rf = 2.3 G and A = 1.34 × 109 ·M for the simulations and Rf = 1.6 G and A = 1.48 × 109 ·M for the experiments. The inset shows the ratio of the resistance of the wide regions of the channel to the resistance of the selectivity filter for the simulations.

It may at first seem counterintuitive that at physiologically relevant ion concentrations the majority of the resistance is in the wide part of the channel. To understand this, consider that the conductance of any medium through which ions flow is a product of the mobility and the charge density. To some extent the effective mobility of the ions is reduced in the selectivity filter by the negative charges of the backbone carbonyl oxygens, yet the concentration is increased enormously (~20 M K⁺), which is more than compensating for the reduction in effective mobility.

	SUMMARY AND CONCLUSIONS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION SUMMARY AND CONCLUSIONS REFERENCES

We have described and carried out a hierarchical strategy for characterizing the permeation of the KcsA potassium channel. This strategy begins with the x-ray crystal structure of the channel (Doyle et al., 1998) from which a series of opened channel structures was built by a process of rotating the transmembrane helices surrounding the channel permeation pathway in a manner generally consistent with blocking experiments, toxin-binding studies, and hydrophobic matching principles. Experimental pH conditions for the functionally closed and open conformations were incorporated through pK_a calculations of the ionization states of the titratable residues. Electrostatics calculations were used to obtain the potential profiles of ions along the channel axis. MD simulations were performed to obtain the ion diffusion coefficient for incorporation into BD calculations for generating ion fluxes. Although the calculations have many approximations (e.g., representing the ion-channel potential profile by an electrostatic calculation, representing the ion motion through the channel as one-dimensional, representing ion-ion interaction in the channel as governed by an equivalent dielectric constant, not having a high-resolution structure for the open form of the channel), the calculations appear to capture the essence of the permeation through the KcsA channel, based on the close correspondence between the behavior of the real channel and the BD simulations.

Examination of the calculations suggests several features of K⁺ permeation through potassium channels. 1) The x-ray model structure of the KcsA channel is closed not by physical occlusion but rather by image forces that arise in attempting to move an ion from the electrolyte to a narrow, nonpolar cavity. 2) Changing the ionization states of the titratable residues by reducing the pH from 7 to 4 in the x-ray crystal structure does not alter the channel-ion potential profile sufficiently to account for the increase in conductance in lowering the pH. Therefore a conformational change in the transmembrane helices of the channel is required for ion conduction. 3) The local friction of ions moving within the channel is similar to that of ions in bulk water. However, the effective mobility in the selectivity filter is smaller due to the partial confinement of ions in the selectivity filter by negative structural charges in the protein (Chiu et al., 1993; Allen et al., 1999). This can be seen in the computed trajectories in Fig. 4, noting that the fluid dynamic diffusion coefficient for ions in the selectivity filter is the same as for ions elsewhere in the channel, but the ions in the selectivity filter are effectively immobilized by the negative charges in the carbonyl oxygens lining the selectivity filter lumen. 4) The most common occupancy state of the open channel at near-physiological concentrations is two ions in the selectivity filter. In and near the selectivity filter are where ion-ion interactions primarily occur. The dominant mode for ion motion through the selectivity filter is effectively by a knock-off mechanism involving three ions, as suggested experimentally (Jiang and MacKinnon, 2000). 5) For near-physiological electrolyte concentrations there is seldom multiple occupancy in the extracellular vestibule or along the intracellular half of the channel. Ion motion in those regions is therefore largely governed by the electrodiffusion of single ions. It is for this reason that the one-dimensional approximation to the motion embodied in the Brownian computational model in this paper gives reasonably accurate estimates of the flux. 6) The overall saturation behavior of the channel in symmetric solutions is rather well fit by a simple model of two resistors in series, one for the selectivity filter and the other for the wider parts of the channel. The resistance of the selectivity filter is ~2 × 10⁹ . The combined resistance of the wider parts of the channel is approximately inversely proportional to the bath concentration of permeant ions.

Our calculations are generally consistent with other computational studies of the KcsA channel. Our MD trajectories for ions in the selectivity filter and in the wider regions of the channel in general agree with those of other workers (Bernèche and Roux, 2000; Guidoni et al., 2000; Shrivastava and Sansom, 2000) The observation in the BD calculations that the dominant occupancy state in the channel is two ions ~7 Å apart agrees with free energy perturbation studies (Åqvist and Luzhkov, 2000), as well as being also suggested by the x-ray model structure (Doyle et al., 1998). Our ionization state for the residue Glu-71 is somewhat more negatively charged than suggested by others (Roux and MacKinnon, 1999; Åqvist and Luzhkov, 2000), perhaps due to a somewhat different methodology of computing the pK_a. As these residues are over 6 Å from the channel axis and are largely occluded by the selectivity filter from the channel axis, this difference is unlikely to affect the conclusions of this study. We also appear to show somewhat less tendency for an ion to occupy the wide part of the channel than previously suggested (Roux and MacKinnon, 1999). This difference may be due to the fact that our ion distributions are computed for an opened channel. Also, Roux and MacKinnon based their main conclusions on calculations in which the protein was assigned a dielectric constant of 2, whereas our assumed dielectric constant of 20 would somewhat attenuate the effect of the pore helix potentials to attract an ion into the cavity.

In the approach described in this paper, molecular modeling based on structural data, electrostatic calculations, MD, and BD simulations are combined to achieve a comprehensive view of the permeation process in the KcsA ion channel. The same overall approach seems capable of application to other ion channels as well.