INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS REFERENCES

With the recent determination of the crystal structure of the KcsA potassium channel (Doyle et al., 1998), the focus of theoretical studies of ion permeation has shifted from gramicidin A to KcsA. So far, the main thrust of these studies has been to employ the atomic structure of the channel in molecular dynamics (MD) simulations to investigate permeation properties such as ion binding sites, selectivity, and diffusion (Guidoni et al., 1999, 2000; Allen et al., 1999, 2000a,b; Bernèche and Roux, 2000; Roux et al., 2000; Åqvist and Luzhkov, 2000; Luzhkov and Åqvist, 2001; Shrivastava and Sansom, 2000; Biggin et al., 2001). Although these MD-based investigations have given useful insights about the function of the KcsA channel, they are limited to the nanosecond time regime and therefore cannot predict the channel currents that require microsecond-long simulations. Because conductance is the main physiological observable characteristic of a channel, it is essential to devise a method that enables its computation. Using a simplified KcsA channel, we have shown in a previous article that three-dimensional Brownian dynamics (BD) simulations can fulfill this task (Chung et al., 1999). The shape of the model channel used in that work resembled that deduced from crystallography, but atoms constituting the channel were represented as a homogeneous, low-dielectric medium with a smooth water-protein boundary. This simplified model was a first step in BD studies of the KcsA channel, and it has provided considerable insight into permeation process in multiple-ion channels. The next step is to incorporate the full structural information obtained from crystallography into BD simulations, so that channel function can be directly related to the structure of the KcsA protein. Here we carry out BD simulations on a potassium channel model that includes all the experimentally determined channel protein, composed of 396 residues and 3504 atoms excluding polar hydrogens.

Paramagnetic spin resonance studies (Perozo et al., 1998, 1999), associated with the pH-dependent gating of this channel (Heginbotham et al., 1999; Cuello et al., 1998), have indicated that the experimental structure corresponds to a closed conduction state and that the transmembrane helices forming the intracellular pore move away from the channel axis during gating. Also, in the crystal structure, the radius of the narrowest section of the pore on the intracellular side is ~1.2 Å, smaller than the radius of the potassium ion (1.33 Å). We therefore first modify the intracellular entrance of the channel with MD simulations to create a set of possible open-state structures. To gain insight into the permeation mechanism, we construct potential energy profiles experienced by an ion traversing the channel that is occupied by other ions. These multi-ion profiles reveal equilibrium ion positions, the effect of charges on ionizable residues required for the channel to conduct, and the rate-limiting step for ions passing the channel. To quantify these predictions, we perform BD simulations, employing various findings from atomic detail MD calculations, namely, ionic diffusion estimates and ionic potential of mean force interactions. By carrying out these simulations we observe the combined effects of gate size and charge distribution on channel conductance, propose a conducting-state channel based on MD simulations with the full channel protein, and improve our understanding of permeation in multiple-ion channels.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS REFERENCES

We employ a combination of fully microscopic MD and semi-microscopic BD simulations to provide connections between the known atomic structure of the KcsA potassium channel and its observable properties. First we describe our methods for generating possible open-state structures with MD and then present our BD technique.

Molecular dynamics: creating an open state

We increase the radius of the intracellular gate in a series of MD simulations, with the CHARMM v.25 package (Brooks et al., 1983), to produce a set of conducting channel candidates to be tested with BD. During MD simulations, the protein is kept neutral (Lazaridis and Karplus, 1999) and the CHARMM 19 extended atom parameter set is employed. The experimental KcsA protein structure (Doyle et al., 1998), taken from structure file 1BL8 of the Protein Data Bank, is modified to include terminal residues and missing side chains, as described previously (Allen et al., 2000a). The structure is then hydrated with SPC/E water molecules and attached to intracellular and extracellular reservoirs of water. Constraints representing the interactions of the protein with a membrane layer are included so that the protein undergoes no large distortions while ion selectivity is unaffected (Allen et al., 2000a). Following 200 ps of equilibration, three potassium ions are placed near the experimentally determined positions, based on Rb⁺ difference maps (Doyle et al., 1998), and held with 100 kT/Ų harmonic constraints (1 kT = 4.11 × 10²¹ J at T = 298 K). The ions are placed inside the wide cavity, near the T₇₅ carbonyl oxygen (a broad maximum in charge density), and near the Y₇₈ carbonyl oxygen. BD simulations, described in Results, show that an ion always resides near the intracellular entrance of the channel due to the presence of glutamic acid residues. Because the presence of a potassium ion will polarize the protein, we also hold an ion in this region during MD simulations.

In BD simulations, the channel protein is assumed to have a rigid structure corresponding to the average positions of atoms forming it. Because the selectivity filter is quite narrow and contains two ions, it is likely to exhibit sizable fluctuations during MD simulations. Therefore, we impose some uniformity on the filter protein by applying a cylindrical harmonic potential (100 kT/Ų) to carbonyl oxygen atoms (residues 75-78) so as to maintain a pore radius of ~1.4 Å. This filter radius is similar to that observed in the crystallographic structure (Doyle et al., 1998), and the overall change to the shape of the filter segment is minimal. Despite this necessary simplification imposed on the model, BD captures the salient conduction properties of potassium channels. This is because the most essential features that govern the permeation of ions across a narrow pore are captured in our model. These are the overall channel shape, the magnitude of charges that creates an energy well deep enough to attract three ions in and in the vicinity of the selectivity filter, the radius of the intracellular gate, and the charges placed on the ionizable residues guarding the entrances. There are small details that can be neglected in the model for the purpose of BD simulations. For example, small variations in the radius of the filter have no perceptible effects on electrostatic calculations and BD results. In contrast, such variations would have a drastic effect on MD results. We also note that the outermost of the two resident ions is ejected from the filter ballistically and almost instantaneously (see, for example, Allen et al., 1999) when the stable equilibrium is disrupted by a third ion penetrating the filter. Thus, the motions of the carbonyl oxygens lining the filter owing to the rapid movement of potassium ions through the conduit will make no difference on the macroscopic properties deduced from BD simulations.

Initially, atomic constraints of 20 kT/Ų are applied to the -carbon atoms in the inner and outer helices. In a series of 100-ps simulations, the target restraints for z < 0 Å are moved away from the channel axis in small steps until the minimum radius of the intracellular pore reaches the desired value. Intracellular gates of approximate minimum radii R_min = 2, 3, 4, and 5 Å are generated with this procedure. All target restraints on the inner and outer helices are then removed and a repulsive cylinder is placed inside the intracellular pore such that any atom whose center resides within z < 0 Å and r < R + 2 Å will experience a harmonic potential (200 kT/Ų). This leaves both the inner and outer helices free to rotate and reach a minimum in the potential surface for this widened structure. It has been suggested that rotation of the transmembrane helices occurs during gating (Perozo et al., 1998, 1999). However, because this evidence is only qualitative, such modifications to the structure have not been imposed during modeling. Also, because the channel is widened by a maximum of ~4 Å, rotation of helices is expected to be small. Simulation with the transmembrane helices released is carried out for an additional 50 ps during which little change in the pore shape occurs because of the placement of the repulsive cylinder. A sample from the trajectory is taken that maximizes the current in BD simulations. Water molecules and ions are discarded and a single subunit is chosen and replicated to impose fourfold symmetry about the channel axis.

Fig. 1 A illustrates the KcsA potassium channel protein structure for an open-state channel with the minimum intracellular gate radius of R_min = 3 Å, hydrated with SPC/E water molecules. For clarity, we show only two of the four subunits that form the pore. The dielectric interface between protein and water is determined by assigning the protein atoms Born radii (Nina et al., 1997) and tracing the channel pore with a water molecule sphere of radius 1.4 Å. The profile obtained from the minimum boundary radius at each axial position is rotated around the channel axis by 360° (Fig. 1 B) to generate a three-dimensional pore. Finally this axially symmetrical pore shape is superimposed on the protein structure. Because of the fourfold symmetry of KcsA, and the large size of the water molecule in comparison to the size of the KcsA pore, the axially symmetric boundary is found to be a good approximation. For the model system illustrated in Fig. 1 B, the boundary has an intracellular pore (23  z < 3 Å) with radius  3 Å, reaching a maximum of 5.5 Å in the cavity region (3  z < 8 Å), whereas that in the selectivity filter (8  z  23 Å) is 1.4 Å. The total length of the channel, including the inner and outer vestibules, is 68 Å. When performing BD simulations, cylindrical reservoirs are attached to the channel ends as shown in Fig. 1 C. Because the reservoirs in our system extend only a short way from the pore, it is the channel protein itself that is forming the dielectric boundary, and not a lipid bilayer. Thus the thickness of the low dielectric slab is larger than that of a typical membrane. The origin (z = 0) is chosen to lie near the center of the wide cavity region. We refer to Doyle et al. (1998) for a guide to the positions of each residue within the KcsA protein structure.

View larger version (59K): [in this window] [in a new window]   FIGURE 1   Model potassium channel. (A) Two of the four subunits of the full experimentally determined protein and the positions of water molecules inside of the pore are illustrated. The aspartate-arginine pairs near the extracellular entrance of the channel (right-hand side) and the glutamate-arginine pair nears near the intracellular entrance of the channel are indicated in red and silver. (B) An outline of the water-protein boundary of a channel with the minimum intra-pore radius of Rmin = 3 Å is shown. The boundary is first smoothed with a water molecule, and then the curve is rotated by 360° to obtain a three-dimensional channel structure. (C) For BD simulations, a reservoir of 30-Å radius is attached at each end of the channel and a fixed number of potassium and chloride ions are placed in each reservoir. The height of the reservoir is adjusted to obtain a desired ionic concentration.

The widening of the intracellular pore by spreading of transmembrane helices is evident in Fig. 2, where pore radius profiles of open-state models with R_min = 3, 4, and 5 Å are compared with that of the crystallographic structure (broken line). We note that these curves are calculated with the use of atomic Born radii of protein atoms (Nina et al., 1997), and not van der Waals radii. As a result the crystallographic structure has a pore radius as small as 0.6 Å in the intracellular gate region. Rotation of the helices around axes parallel to the z axis during modeling of the R_min = 3 Å channel never exceed 5°. The axes chosen for this measurement, for the inner and outer helices, pass through the C₉₀ and A₅₀ -carbon atoms, respectively.

View larger version (19K): [in this window] [in a new window]   FIGURE 2   Pore radius profiles of four different open-state channels and that of the crystallographic structure (  ). The numbers accompanying the traces are the minimum pore radii (Rmin) of the entrance from the intracellular side. Pore radii of Rmin = 3, 4, and 5 Å are shown, whereas the channel with Rmin = 2 Å has been omitted for clarity.

Brownian dynamics simulations

We perform three-dimensional BD simulations to predict the channel conductance. The method of implementing BD simulations is detailed elsewhere (Li et al., 1998; Chung et al., 1998, 1999). Here we give a brief description of the method and refer to our previous publications for further details. In BD simulations, the trajectories of ions are followed by integrating the Langevin equation:

(1)

where m_i, v_i, m_ii and q_i are the mass, velocity, friction coefficient, and charge of an ion with index i. The quantities F_Ri, E_i, and F_Si are the random force, systematic electric field, and short-range force experienced by the ion, respectively. Electrostatic forces are determined by solving Poisson's equation with the boundary element method (Hoyles et al., 1996, 1998b) and implemented in BD simulations using look-up table techniques (Hoyles et al., 1998a).

The Langevin equation (Eq. 1) is solved using the algorithm of van Gunsteren and Berendsen (1982). A multiple time-step algorithm is used, where a time-step of t = 100 fs is employed in the reservoirs and 2 fs in the channel. The latter is necessitated by the fact that forces change rapidly inside the channel. Recalling that the relaxation time constant, ¹, is 31 ps for potassium ions, their motion corresponds to over-damped diffusion in the reservoirs (t ¹) but not in the channel. A temperature of 298 K is assumed throughout. The friction coefficient m_ii in Eq. 1 is related to the diffusion coefficient D_i by the Einstein relation D_i = kT/m_ii. Bulk ionic diffusion coefficients (1.96 × 10⁹ m²/s and 2.03 × 10⁹ m²/s for K and Cl ions, respectively (Lide, 1994)) are employed in reservoirs, but reduced values are used in the channel according to the MD results (Allen et al., 2000b). Within the intracellular pore, we maintain a bulk diffusion coefficient for the K ion, whereas a value of 50% bulk diffusion is employed in the hydrophobic cavity region. Because the value of the diffusion coefficient used in the selectivity filter has been shown to have little effect on channel conductance (Chung et al., 1999), we also use a value of 50% bulk in that region.

The dielectric constant of the protein is assumed to be p = 2 in the main body of this work, though we do consider the effect of small variations (2-5) on the results. The situation with regard to the channel water is less certain. It is expected that the dielectric constant of the channel water will be reduced from its bulk value of 80 (Sansom et al., 1997), but the amount of reduction is highly dependent on the process one uses. For example, the response of water molecules to an external perturbation is highly suppressed, so this type of calculation would result in a very low  value. On the other hand, if one considers dielectric screening of an ion, a much larger  value is expected because the potassium ions appear to be well solvated throughout the channel (Allen et al., 1999, 2000a). In past simulations on a model KcsA channel (Chung et al., 1999), a value of 60 was found to optimize conductance, which we adopt here as well. Because the boundary element method does not permit ions to cross dielectric boundaries, a uniform dielectric constant of 60 is employed throughout channel and bath water. The change from reservoir (80) to channel water (60) is represented by a Born energy barrier of

(2)

which assumes a Born radius of R_B = 1.93 Å, determined from the enthalpy of hydration for K⁺ ions (Bockris and Reddy, 1970). This barrier is introduced as a smooth switching function erected at the channel entrances (Chung et al., 1999). The short time-step of 2 fs used in the channel pore (25  z  25 Å) ensures that ions are subjected to these short-ranged forces. Ions also interact with the channel protein lining via a repulsive short-range potential similar to that described in Chung et al. (1999). The Coulomb interaction between ions is modified at short distances by strongly repulsive contact, weakly attractive dispersion, and oscillating hydration forces. A correct modeling of these short-range contributions is very important in multi-ion channels because they have a direct bearing on the relative positions of ions within the channel. This has been highlighted in recent MD simulations of the KcsA channel that revealed the role played by hydration waters in positioning of ions in the selectivity filter (see, for example, Åqvist and Luzhkov, 2000; Allen et al., 2000a; Bernèche and Roux, 2000). We have previously incorporated these short-range ion-ion interactions in BD simulations by parameterizing the potentials of mean force obtained from MD simulations (Corry et al., 2001). We adopt the same analytical form here:

(3)

where U₀ is the strength of the potential, R_c is the contact distance defining the repulsive first term. The second term describes the damped oscillations due to hydration forces with R  R_c, c_e = 1 Å, and c_w = 2.76 Å (water diameter). Table 1 lists the parameters for each ion pair.

A cylindrical reservoir of 30-Å radius and variable length is connected to each end of the channel (see Fig. 1 C). We note here that the reservoir boundaries simply serve to confine the ions within the simulation system, thus maintaining the average concentrations in the baths at the desired values. Scattering of an ion from the boundary wall can be viewed as an ion moving out of the reservoir and another one entering at the same time. Implicit in the use of reservoirs is the assumption that electrolyte solution continues beyond its boundaries. Symmetric concentrations are generated with 16 potassium and 16 chloride ions in each reservoir, and the reservoir length is adjusted to set the concentration. With the exception of the current-concentration curves, 300 mM is used in all simulations, which corresponds to a reservoir length of 32 Å. During conduction events, the average concentration in the reservoirs is maintained with a stochastic boundary (Chung et al., 1999).

In experiments, a potential difference is applied via electrodes at macroscopic distances from the channel. In a small system such as the one used in BD simulations, this procedure cannot be simply mimicked by fixing the potential at the boundaries. Because the channel is highly charged, this would distort the potential in the channel vicinity compared with the experimental set-up. To circumvent this problem, we generate the potential difference across the channel via a uniform applied electric field of strength E₀. Operationally, this is equivalent to having a pair of isopotential plates located at large distances from the channel. In the absence of any dielectric boundary, the potential difference across a channel of length l would be simply lE₀. The presence of a dielectric boundary and charges on the protein wall, however, distorts this field in and nearby the channel. Thus, the precise potential difference across the channel depends on the selected reference points at the two sides of the potassium channel. A convenient choice for our purposes is the center of each reservoir at z = ±50 Å, which we employ throughout this work. To determine the potential difference across the reference points, we measure sample potentials along the z axis once every 10 ps and average over the simulation period with 300 mM KCl in the reservoirs. The difference in the average potentials between the reference points, V, is taken to be the potential difference across the channel. In the following, unless otherwise stated, we apply an electric field of strength E₀ = 2 × 10⁷ V/m in BD simulations. The corresponding values of V for each channel are given in Table 2. Note that because of the highly asymmetric charge distribution in the KcsA channel, the V values are not symmetric.

The energy profile of an ion along the z axis is constructed by solving Poisson's equation with a test ion fixed at regular positions along the channel axis. The energy profile seen by an ion entering the channel already housing one or more cations is found using a steepest-descent algorithm (Chung et al., 1999), where forces on all resident ions are minimized and the total energy of the system is recorded at each test ion position.

The BD program is written in FORTRAN, vectorized, and executed on a supercomputer (Fujitsu VPP-300). The amount of vectorization varies from 67% to 92%, depending on the number of ions in the reservoirs. With 64 ions in the reservoirs, the CPU time needed to complete a simulation period of 1 µs is ~53 h. Simulations under various conditions, each lasting 1 × 10⁶ time steps (0.1 µs), are repeated 10-30 times.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUDING REMARKS REFERENCES

Charges on ionizable residues

The KcsA sequence (Doyle et al., 1998) consists of a number of ionizable residues close to the permeation pathway that may affect channel conduction. These residues, and their approximate positions, are listed in Table 3. Our model channel is taken to be overall neutral by grouping ionizable residues into likely dipoles. The presence of any net local charge severely hampers ion permeation, presumably due to the low dielectric strength of the protein. As shown in Table 3, there is only one unpaired arginine residue R₂₇ in the protein sequence. This residue is taken to be neutral in our model because it is found to have a strong destabilizing influence on permeating ions in this sector of the channel (Allen et al., 2000a). We stress that R₂₇ is in a region where 21 residues from the N-terminus and 39 residues from the C-terminus are missing from the experimental KcsA sequence, and also it is in close proximity to the lipid bilayer surface. Therefore, it is reasonable to assume that R₂₇ could be neutralized by missing charged residues or ions in solution.

The E₅₁/R₅₂ pair is far from the channel axis and has been found to have little impact on conduction. These residues are therefore kept neutral. For reasons to be discussed below, the E₇₁/R₈₉ pair are also kept uncharged in our simulations. The remaining ionizable residues D₈₀/R₆₄ and E₁₁₈/R₁₁₇ are required to be charged for the channel to conduct. We have examined the influence of these dipoles on the conductance properties of the R_min = 3 Å channel by performing BD simulations. We assign charges q_e and +q_e on the D₈₀/R₆₄ pair near the extracellular side, creating four mouth dipoles. Similarly, charges q_i and +q_i are placed on E₁₁₈/R₁₁₇ residues to create four intracellular mouth dipoles. These charges are distributed among the atoms of ionizable side chains, following the method described by Lazaridis and Karplus (1999), and have been indicated in Fig. 1, A and B. Both set of dipoles provide attractive energy wells for cations to enter the pore and prevent anions from entering it. Analysis of MD trajectories reveals that the fourfold symmetry of the E₁₁₈ side chains is maintained in the presence of an ion in the vicinity of these residues. Thus the rigid tetramer employed in BD simulations is appropriate and does not introduce an artificially deep energy well at the intracellular entrance.

We show in Fig. 3 the variation in outward (filled circles) and inward (open circles) channel currents with charge q_i on intracellular dipoles (Fig. 3 A) and q_e on extracellular dipoles (Fig. 3 B). The positions of these polar residues are indicated in the inset as dark and white balls. Electric fields of +E₀ and E₀ are applied to generate the outward and inward currents, respectively (see Table 2 for the corresponding potential differences). The charge q_e is fixed at 0.9e during optimization of q_i, whereas q_i is fixed at 0.7 e during optimization of q_e. When q_i is instead held at 0.3 e during optimization of q_e, very similar outward flowing currents are observed, as shown in Fig. 3 B (gray triangles). Both outward and inward currents increase rapidly with q_i and remain high between q_i = 0.3 e and 0.8 e (Fig. 3 A). The current declines slightly when the charge is further increased to 0.9 e. In contrast, both outward and inward currents reach a maximum when q_e is near the full unit charge e (Fig. 3 B). A small deviation in q_e is seen to cause a rapid decline in the channel current.

View larger version (23K): [in this window] [in a new window]   FIGURE 3   Dependence of channel current on fixed charge strengths for channel with Rmin = 3 Å. (A) The current passing through the channel under a driving field of E0 () or E0 () is plotted against the charge qi on each of the four glutamate (E118) and arginine (R117) pairs. The charge on each of the four aspartate (D80) and arginine (R64) pairs are held constant at 0.9 e during the optimization of qi. (B) Same as in A but against the charge qe on each of the four aspartate (D80) and arginine (R64) pairs. The charge qi is fixed at 0.7 e during optimization of qe. Similar results are obtained when the charge qi is fixed at 0.3 e during optimization of qe, with the outward current shown as gray triangles. In this and all subsequent figures, a data point represents the average of 10-20 sets of simulations, each set lasting 1 × 106 time steps (or 0.1 µs). Error bars in this and following figures have a length of 1 SEM and are not shown when they are smaller than the data points.

Fig. 4 shows the variation of outward (filled circles) and inward (open circles) currents with charge q_i on the intracellular dipoles for the channels with R_min = 2, 4, and 5 Å (A, B, and C). The relative magnitude of the current is seen to be sensitive to q_i for R_min = 2 Å, but this dependence is much less pronounced for the larger channels with R_min = 4 or 5 Å. The curves for optimization of q_e show similar trend as those illustrated in Fig. 3 C, suggesting use of full unit charges for optimal currents.

View larger version (14K): [in this window] [in a new window]   FIGURE 4   Dependence of channel current on fixed charge strengths and channel radius. The charge on the extracellular dipoles (D80/R64) are fixed at qe = e throughout. (A) The outward () and inward () currents across the channel with Rmin = 2 Å under a driving field of ±E0, are plotted against the charge qi on the intracellular dipoles (E118/R117). (B) Same as in A but for a channel with Rmin = 4 Å. (C) Same as in A but for a channel with Rmin = 5 Å.

In the preceding series of simulations, the ionizable residues E₇₁ and R₈₉ are kept neutral. Despite the proximity of the E₇₁ side chain to the selectivity filter, calculations have shown that it is energetically favorable for E₇₁ to be protonated (Roux et al., 2000; Luzhkov and Åqvist, 2000). However, recent calculations have predicted that the presence of multiple cations can increase the probability of this residue being charged (Ranatunga et al., 2001). Here we examine the effect of placing charges on the E₇₁ and R₈₉ residues. First, we shift a fraction of the unit charge e on the D₈₀ residue to the E₇₁ residue (i.e., the proton on E₇₁ spends a fraction of its time on D₈₀) and measure the current across the channel. This process is repeated until the full charge is transferred from one residue to the other. In this first series of simulations, the R₆₄ residue is fully protonated, whereas residue R₈₉ is in its neutral state. In Fig. 5, the current is plotted against the charge on the E₇₁ residue (filled circles). When the full charge e is on the D₈₀ residue, the outward current under the driving field of +E₀ is at a maximum. The current rapidly decreases as the charge is shifted from D₈₀ to the E₇₁ residue, reaching a minimum at halfway through and remaining at this 10% level thereafter. We repeat the same series of simulations with the full charge on R₈₉ and no charge on the R₆₄ residue. The outward current across the channel first increases slightly and then declines rapidly as the charge on the aspartate residue is gradually transferred from the D₈₀ residue to the E₇₁ residue (open circles in Fig. 5). With the positive charge shifted to R₈₉, the largest number of potassium ions are translocated when 0.2 e is placed on E₇₁ and the remaining 0.8 e on the D₈₀ residue. The current rapidly attenuates, as before, when the charge on the glutamate residue is further increased. We choose R₆₄ protonated and R₈₉ deprotonated to obtain the maximum current, as seen in Fig. 5. With both R₆₄ and R₈₉ protonated and both D₈₀ and E₇₁ carrying the full charge e, the energy well in the pore becomes deep enough to attract on average 4.5 ions in the absence of an applied electric field. Although the system could be made to conduct, this observed number of resident ions in the filter and wide cavity is inconsistent with the crystallographic observations (Doyle et al., 1998). We thus assume that one of the arginine residues, possibly R₈₉, is locally neutralized by a foreign anionic species and only one charge is shared between D₈₀ and E₇₁, but predominantly on the former in the conducting state.

View larger version (12K): [in this window] [in a new window]   FIGURE 5   Variation of currents with the relative charge placed on the two ionizable residues (D80 and E71) near the extracellular channel entrance for Rmin = 3 Å. A fraction of the charge on D80 is transferred to E71, and the current across the channel is determined. The total charge shared between the two residues is kept constant at e. In one series of determinations, the full positive charge is placed on the R64 residue, while keeping the R89 residue neutral (). The same series of simulations is carried out once again, placing the full positive charge on the R89 residue and keeping the R64 residue neutral ().

We conclude from these results that a value of 0.9  q_e  1.0 e for the D₈₀ and R₆₄ residues, and 0.3  q_i  0.8 e for the E₁₁₈ and R₁₁₇ residues provide near optimal channel performance in both directions of permeation for R_min = 3 Å. Also, the E₇₁ residue near the selectivity filter must be kept neutral for conduction to take place. In the following simulations, we consider two combinations of mouth dipoles with q_i = 0.3 and 0.7 e (the former value is assumed if not explicitly stated). For the 2-, 4-, and 5-Å channels, q_i = 0.5 e is used. The full unit charges are placed on the extracellular dipoles (D₈₀ and R₆₄) throughout. We remark that the KcsA channel is usually closed at neutral pH but can remain in a conducting state at low intracellular pH (Heginbotham et al., 1999). Because the E₁₁₈ side chains are in close proximity to the cytoplasmic side of the membrane and accessible to channel water, they are subject to protonation at low pH. We conjecture that the open state of the KcsA channel is likely to involve a reduced charge on the intracellular mouth dipole. Had we placed the full electronic charges on the glutamate and arginine residues guarding the entrance of the intracellular pore throughout, the outward current would have been reduced on average by 46% and the inward current by 26% from the values we report here.

Energy profiles

We illustrate in Fig. 6 A the potential energy of a potassium ion brought into the channel from the intracellular reservoir in the absence of other ions. The charge on the intracellular dipole is fixed at q_i = 0.7 e. Profiles with q_i = 0.3 e (not shown here) are similar to those illustrated in Fig. 6, except that the well near the intracellular entrance of the channel is shallower. The energy well experienced by a single ion is deep, reaching 67 kT within the selectivity filter with respect to the intracellular reservoir. Fig. 6 B shows the energy profile experienced by a second ion entering the channel already housing one ion. The position of the resident ion when the test ion reaches the center of the channel (z = 0 Å) is indicated as an upward arrow. In constructing this and subsequent profiles, we allow the resident ions in the channel to adjust their positions so that both axial and radial components of the force on them are minimized. Thus, the potential profile for multiple ion systems, calculated with a steepest-descent algorithm (Chung et al., 1999), corresponds to the total electrostatic energy required to bring in the charge on the test ion from an infinite distance in infinitesimal amounts. With two ions in the selectivity filter, a third ion entering the channel sees an energy well, created by the mouth dipoles near the intracellular entrance (Fig. 6 C). The ion needs to overcome a central barrier of 7-8 kT to move toward the wide cavity region. Once in the cavity, it experiences only a shallow energy well signifying the low stability of the three-ion system. A fourth ion entering the channel sees a small energy well and then a steeply rising energy barrier within the intracellular pore (Fig. 6 D). The three-ion equilibrium is disrupted when the fourth ion enters the pore, resulting in expelling of the outermost ion.

View larger version (10K): [in this window] [in a new window]   FIGURE 6   Electrostatic energy profile of an ion traversing the channel with Rmin = 3 Å. The potential energy of a potassium ion is plotted in 1-Å intervals in the z direction with a driving field of E0/2 (V = 118 mV). The ion's z coordinate is fixed, but it is allowed to move to its minimum energy position in the x-y plane. In the absence of any resident ions in the channel, the test ion encounters a deep energy well of a depth of 67 kT at z = 13 Å (A). The energy profiles are constructed in the presence of 1 (B), 2 (C), and 3 (D) ions in the channel. The upward arrows indicate the location of the resident ions when the test ion is at the center of the channel (z = 0 Å). The schematic channel in the inset shows the positions of the ions in case C.

The electrostatic profiles constructed above ignore the influence of ionic atmosphere and other dynamic effects such as entropy. BD simulations indicate that anions do not enter the channel. Therefore, neglect of the counterions is justified, and one can still gain useful insights on the permeation process by examining the potential profiles in the channel. From the profiles illustrated in Fig. 6, it is clear that conduction will not take place unless there are two or three resident ions in the channel. For a single ion to traverse the selectivity filter, it has to gain an energy of nearly 60 kT to exit the energy well shown in Fig. 6 A. Similarly, the energy well seen by a second ion is ~40 kT (Fig. 6 B), also too deep to allow permeation. In the presence of two resident ions (Fig. 6 C), the third ion needs to overcome a smaller energy barrier of ~7-8 kT, when moving from the intracellular entrance up into the cavity. This is expected to be the rate-determining step for conduction. These suppositions are tested in the following sections using BD simulations.

Electric potential across the pore

It is of interest to determine the electric potential profile inside the channel during BD simulations. To construct such a profile, we measure once every 10 ps the electric potential along the z axis of the channel in 1-Å steps, and average over the simulation period. The average potential profiles created by applied electric fields of +E₀ (solid line) and E₀ are depicted as solid curves in Fig. 7, A and B, respectively, for a channel with minimum gate radius 4 Å. The uniform electric field, which may be imagined to be generated by a pair of isopotential plates located at large distances, is distorted by the dielectric boundary. The potential differences generated between z = ±50 Å by these two applied electric fields are listed in Table 2. As mentioned in Methods, we take the measured potentials in plotting the current-voltage curves (see later). Inside the pore, the potential does not change linearly; instead, it fluctuates wildly owing to the presence of the fixed charges on the protein wall and resident potassium ions in the channel. The prominent oscillations seen in the selectivity filter are due to the two resident ions, whose locations are shown in Fig. 8 B. Superimposed on the curves are the potential profiles obtained in the absence of ions and protein atom charges (broken lines). These profiles, which are similar to those found in Poisson-Boltzmann solutions (Roux et al., 2000), illustrate that the potential in the pore changes rapidly in the narrow segment of the channel and more slowly in the wider segment.

View larger version (12K): [in this window] [in a new window]   FIGURE 7   The electric potential profiles in the channel with Rmin = 4 Å. The potential along the z axis of the channel is averaged over 10-ns BD simulations and drawn as solid curves. All atomic partial charges in the protein are included, and no restrictions are placed on ionic positions. Uniform electric fields of E0 (A) and E0 (B) are applied across the channel to create the potential differences between the centers of the reservoirs (see Table 2). The potential profiles in the absence of protein atomic charges and ions are drawn as dashed curves for comparison.

View larger version (14K): [in this window] [in a new window]   FIGURE 8   Average number of ions in the channel with Rmin = 3 Å with no applied field (A) and with an applied field of E0/2 (B). The channel with qi = 0.7e is divided into 100 sections, and the average number of ions in each section is calculated over a simulation period (0.1 µs). The outline of the channel and the approximate locations of ions in the absence of an applied field are shown in the inset.

Ions in the channel

In Fig. 8, the histograms of ions in the channel with R_min = 3 Å are illustrated. To construct the histograms, the channel is divided into 100 thin sections, and the average number of ions over a 0.1-µs simulation is plotted. The concentration of KCl in the reservoirs are kept constant at 300 mM and q_i = 0.3 e. With no applied potential, there are three regions in the selectivity filter and cavity where ions dwell preferentially (Fig. 8 A). The locations of their maxima are z = 4.5 Å (inside the cavity), 12.3 Å (near the T₇₅ carbonyl oxygen), and 19.2 Å (near the Y₇₈ carbonyl oxygen). Also, there is another prominent peak in the histogram, centered at 21 Å, near the E₁₁₈ side chain. The average number of ions is 2.9 in the selectivity filter and the cavity and 0.9 near the intracellular entrance of the channel. The ion in the cavity region tends to reside mostly off-axis, as a result of the strong interaction with the oxygen atoms of the T₇₄ residue at the base of the pore helix. These positions are in close agreement with positions observed in Rb⁺ x-ray diffraction maps (Doyle et al., 1998). When an electric field E₀ is applied across the channel so that ions steadily move from inside to outside, the pattern of the histogram shifts, as illustrated in Fig. 8 B. There are now two prominent peaks in the selectivity filter region and one peak, as before, in the intracellular entrance of the channel. The average number of ions in the selectivity filter and cavity is now 2.2, whereas that near the intracellular cellular entrance is 1.1. We remark that with charge q_i = 0.3 e on the intracellular dipoles, the distributions are very similar, except that 1.0 ion resides near the intracellular entrance. In Fig. 9, we illustrate three histograms obtained from the channels with R_min = 2, 4, and 5 Å with no applied field. All histograms show three distinct peaks in the selectivity filter and the cavity. The average number of ions in the channel is similar in each case, 2.7 in the extracellular side and 1.3 in the intracellular side.

View larger version (15K): [in this window] [in a new window]   FIGURE 9   Average number of ions with no applied potential for the channels with Rmin = 2 Å (A), 4 Å (B), and 5 Å (C). For all histograms, qi = 0.5e is used.

Current-voltage relationships

We study the conductance properties of potassium ions under various conditions by performing BD simulations. The current-voltage curves for the channels with different intra-pore radii are shown in Figs. 10 and 11. For the 3-Å channel, two current-voltage curves are presented, one with q_i = 0.7 e (Fig. 10 A) and the other with q_i = 0.3 e (Fig. 10 B). Both curves are fairly linear through the origin when the applied potential is less than ±150 mV, but they deviate systematically from Ohm's law with a further increase in the membrane potential. This nonlinearity results from the presence of an energy barrier in the channel. Intuitively, a barrier is less of an impediment to a permeating ion when the driving force is large. The outward conductance of the channel at an applied potential of 150 mV is 34 ± 4 pS with q_i = 0.3 e, and a similar value is obtained for q_i = 0.7 e. The inward and outward currents are approximately symmetrical when q_i is 0.7 e (Fig. 10 A), but the curve becomes pronouncedly asymmetrical when q_i is reduced to 0.3 e (Fig. 10 B).

View larger version (8K): [in this window] [in a new window]   FIGURE 10   Current-voltage relationships for the channel with Rmin = 3 Å. The magnitude of the current passing through the channel with a symmetric solution of 300 mM KCl in both reservoirs is plotted against the applied potential. The value of qi placed on the E118/R117 residues is 0.7 e (A) and 0.3 e (B).

View larger version (8K): [in this window] [in a new window]   FIGURE 11   Current-voltage relationships for the channels with Rmin = 2 Å (A), 4 Å (B), and 5 Å (C). The simulation conditions used to generate the data points are the same as in Fig. 10.

The current-voltage curves obtained from the channels with R_min = 2, 4, and 5 Å are illustrated in Fig. 11. The features observed in the 3-Å channel are also apparent in the curves shown in Fig. 10. In each of the three curves, there is a systematic deviation from linearity at high applied potentials, with the degree of nonlinearity being most pronounced in the 2-Å channel. Also, slight asymmetries between the inward and outward currents can be discerned from the curves. The conductances at 150 mV for the channels with R_min = 2, 4, and 5 Å are, respectively, 6 ± 2, 147 ± 7, and 172 ± 15 pS. The corresponding values at 150 mV are 7 ± 1, 96 ± 4, and 93 ± 12 pS.

The current-voltage relationships from the KcsA potassium channel have been determined by several groups. The values of the outward conductance at 100 mV they report are 135 pS with 200 mM KCl (Cuello et al., 1998), 97 pS with 250 mM KCl (Meuser et al., 1999), and 83 pS with 100 mM KCl (Heginbotham et al., 1999). Taking concentration dependence into account (see Fig. 16), comparison of these values with the calculated conductances suggest an intra-pore radius of R_min = 3-4 Å. Note that there are substantial variations in the measured values, and more consistent measurements of the KcsA channel conductance would be desirable to determine the model parameters with more certainty. There are also inconsistencies in the degree of rectification observed in experiments. For example, results of Heginbotham et al. (1999) exhibit a pronounced outward rectification, whereas those of Meuser et al. (1999) show symmetrical inward and outward currents. As both experiments were carried out at similar pH levels, it is difficult to understand the origin of this discrepancy. Our simulations (Fig. 10) suggest that rectification arises from a reduction of charges on the glutamate residues E₁₁₈ as a result of protonation at lower pH levels. Thus, accurate measurement of pH dependence of inward and outward currents would be very valuable in pinpointing the charge state of the E₁₁₈ residues.

Dielectric constant of protein

In the preceding analyses, we used _p = 2 for the dielectric constant of protein. Unlike water and lipid, which form homogeneous media, proteins are quite heterogeneous, exhibiting large variations in polarizability depending on location (interior or surface) and nature of the process (for recent reviews, see Nakamura, 1996; Schutz and Warshel, 2001). From a microscopic point of view, this should make assigning a fixed _p value to an entire protein in continuum electrostatic calculations problematic, but in practice, it seems to work quite well. In the following series of simulations, we use _p = 3.5 and 5 in solving Poisson's equation and compare the results of BD simulations obtained by using these two different values to those obtained by assuming _p = 2. We show that the precise value of the dielectric constant of protein we adopt has only negligible effects on the macroscopic properties one derives from BD simulations.

We first compare in Fig. 12 A the energy profiles encountered by an ion traversing the channel when _p is assumed to be 2 (broken line) and 5 (solid line). In this and the following figures, we use the channel with the R_min = 4 Å. The profiles are obtained with two ions placed in the selectivity filter. With 0.5 e placed on the E₁₁₈ and R₁₁₇ and an applied field of +E₀, the depth of the barrier height V_B encountered by the third ion attempting to proceed toward the inner cavity increases from 2.75 kT to 3.37 kT as _p is changed from 2 to 5. In Fig. 12 B, we show the variation of outward currents under the influence of an applied field of +E₀ with the charge placed on the E₁₁₈ and R₁₁₇ residues (filled circles). With _p = 5, the peak current is observed when the ionizable residues carry 0.6 e. Superimposed on the graph are the data for outward currents illustrated in Fig. 4 B, where the simulations are carried out with _p = 2 (open circles). The current is slightly attenuated when _p is increased, and the optimum charge for the glutamate-arginine pairs shifts from 0.5 e to 0.6 e.

View larger version (13K): [in this window] [in a new window]   FIGURE 12   Energy profiles and the variation of currents with charges on the glutamate-arginine residues obtained with p = 2 and 5. The field of 2 × 107 V/m is applied. The 4-Å-channel model is used for this and the subsequent two figures (Figs. 13 and 14). (A) The profiles are constructed in the presence of two ions in the selectivity filter with p = 5 () and p = 2 (  ). (B) The outward current plotted against qi with p = 5 () is compared with that obtained with p = 2 (), reproduced from Fig. 4 B.

The location and the number of ions in the channel in a conducting state also remain unchanged when _p is increased from 2 to 5, as shown in Fig. 13. The histograms are tabulated during the simulation period of 0.1 µs in the presence of an applied field of +E₀. The average number of ions in the selectivity filter and inner chamber is 2.23 when _p = 5 and 2.33 when _p = 2. The corresponding numbers in the intracellular half of the channel are 1.65 and 1.49. 

View larger version (20K): [in this window] [in a new window]   FIGURE 13   Average number of ions in the channel with the assumed values of p = 2 (A) and p = 5 (B). The histograms are constructed in the presence of an applied field of 2 × 107 V/m.

Fig. 14 shows the current-voltage relationships derived by using _p = 3.5 and 5.0. Superimposed on these curves is the curve reproduced from Fig. 11 B with _p = 2. The currents at all driving voltages are slightly attenuated when the _p is increased. This decrease is due to a small increase in the height of the energy barrier V_B, as illustrated in Fig. 12 A.

View larger version (10K): [in this window] [in a new window]   FIGURE 14   Current-voltage relationships derived with p = 3.5 (A) and p = 5 (B). For comparison, the curve fitted for the current-voltage date illustrated in Fig. 11 B is superimposed in A and B (  ).

There are several microscopic investigations of the dielectric constant of proteins from MD simulations (Smith et al., 1993; Simonson and Perahia, 1995; Simonson and Brooks, 1996; Simonson, 1998; Pitera et al., 2001). The consensus emerging from these MD studies of globular proteins is that the effective dielectric constant for the whole protein varies between 10 and 40 depending on the choice. However, when only the interior region of the protein consisting of the backbone and uncharged residues is considered, _p values drop to 2-4. Thus, the relatively large effective _p values extracted from MD are mostly due the charged side chains on the protein surface. These studies suggest using a low _p value in the protein interior and a transition region between the protein and solvent with a medium  value that takes into account the polarizable side chains as well as the more ordered layer of water surrounding the protein. Estimating _p for proteins forming the channel wall using MD calculations presents an inverse problem. Here, water is embedded in protein rather than protein embedded in water. Although the quantitative results will certainly be different from globular proteins, we expect the general features mentioned above to apply to the membrane proteins. This issue clearly deserves further MD studies. Our studies, nevertheless, indicate that the calculations carried out using _p = 2 remain largely valid even if the effective dielectric constant of channel protein turns out to be substantially higher than 2.

Analyses of conduction processes

In analogy with the flow of a liquid through a constricted pore, one would naively expect that the passage of ions across the narrow selectivity filter is the rate-limiting step for conduction across the potassium channel. Our detailed analysis of conduction process reveals that this is not the case. This is because ions diffuse singly through the intracellular gate and cavity regions but undergo a rapid multi-ion shuttling process within the selectivity filter. We have shown earlier that ionic mobility in the filter region has little impact on conductance because permeation through the filter is governed by the electrostatic repulsion of ions (Chung et al., 1999). These rapid movements of ions in the filter have also been observed in MD studies (Allen et al., 1999); thus, they are not artifacts of using a rigid filter structure in the BD simulations.

We first discuss the potential energy profiles for three K⁺ ions in the channel as they provide an intuitive picture for the permeation process. The profiles shown in Fig. 15 A are constructed as in Fig. 6 C: two ions are placed in the selectivity filter and their positions are varied to minimize their energy as a third ion is moved from z = 50 to z = 10 Å in 1-Å steps. A driving field of E₀ is applied to move ions outward. It is seen that for permeation to take place, an ion in the energy well near z = 20 Å has to go over the barrier at z = 8 Å. The height of this barrier V_B, measured as the potential energy difference between z = 20 and 8 Å, decreases as the intra-pore radius is increased: V_B = 11.1, 6.9, 3.3, and 1.3 kT for R_min = 2, 3, 4, and 5 Å, respectively (here q_i = 0.3 e is used for the 3-Å channel). Note that in all cases, once an ion is over the barrier, it is driven toward the selectivity filter under a potential gradient. The Coulomb force exerted by this ion approaching the cavity disrupts the equilibrium of the two resident ions in the filter (see Fig. 8 B), resulting in ejection of the outermost ion to the reservoir and thus completing the conduction event.

View larger version (13K): [in this window] [in a new window]   FIGURE 15   Rate-limiting steps for ion permeation. (A) The energy profiles presented to a potassium ion as it proceeds toward the cavity from inside the cell are plotted. A schematic diagram in the inset shows the location of two resident ions in the selectivity filter when the permeating ion is at z = 0 Å. The numbers accompanying the profiles indicate the intra-pore radius of the channel. The ion traversing the 2-Å channel, for example, must overcome the barrier whose height is indicated as VB. (B) The potential profiles presented to the innermost ion in the selectivity filter as it moves toward the intracellular space are plotted. The three profiles correspond to channels with Rmin = 2, 3, and 4 Å. The barrier height for the 2-Å channel is indicated as VB. In the inset, the positions of the resident ions and the test ion at its starting point are indicated. The profiles obtained from the 5-Å channel, not shown here for clarity, are similar to that obtained from the 4-Å channel.

The profiles encountered by an ion transiting in the inward direction under a driving field of E₀ are similarly constructed and illustrated in Fig. 15 B. Here, two ions are placed in the selectivity filter and one ion at the bottom of the energy well near the intracellular entrance. A test ion is then moved from z = 10 Å to z = 50 Å, while allowing the other three ions to adjust their positions. The barrier heights encountered by the test ion between z = 8 and 8 Å are given by V_B = 7.7, 4.9, 3.8, and 3.8 kT for R_min = 2, 3, 4, and 5 Å, respectively. The differences between the inward and outward barriers provide a qualitative explanation for the rectification of the I-V curves observed in Figs. 10 and 11. For example, the substantial outward rectification seen in the 5-Å channel (Fig. 11 C) arises from the outward barrier being three times smaller than the inward one.

From these profiles, we expect that an ion attempting to climb over the barrier to move more slowly and spend more time in the deeper wells compared with shallower ones, thus giving a simple explanation for the rapid increase in conductance with radius seen in Figs. 10 and 11. We quantify the above statement through an analysis of the ion trajectories obtained from BD simulations. For reference, we note that the drift velocity of potassium ions in a bulk solution under a driving field of E₀ = 2 × 10⁷ V/m is 1.5 m/s. The average drift velocity of ions in the channel as they climb up the energy barrier in Fig. 15 is found to be substantially smaller than the drift velocity in a bulk solution, but it becomes larger than the bulk value once they are over the barrier and move toward the cavity. For example, for the 2-Å channel, the respective velocities before and after the barrier are 0.24 and 5.73 m/s under the driving field of E₀. Thus, once an ion climbs over the energy well, the transit across the remaining part of the pore is rapid. As the channel gets larger, the barrier becomes less steep and the drift velocities on either side of the barrier gradually approach toward the bulk value. Here, velocities are determined by averaging the ionic speeds during conduction processes whenever there is an ion in the desired segment.

The time an ion spends in various parts of the channel also provides valuable insights into the permeation process. The maximum number of ions the channel can process is expected to be determined primarily by the time it takes for an ion to climb out of the well located near z = 20 Å. In Table 4, we tabulate the average time the well remains empty (t₀) and occupied (t₁) by an ion for the four channels with different intra-pore radii. The mean transit time is essentially the sum of the times t₀ and t₁. The waiting times shown in Table 4 correlate well with the barrier heights (Fig. 15), clearly demonstrating that the rate-limiting step for conduction is passing of the ions over the central barrier. For the most part, an ion resident in the well makes repeated attempts to move toward the cavity, but success is rare and is proportional to the barrier height. Analysis of BD trajectories reveals that, under the applied field of E₀, the success rates are 1.5%, 4%, and 5.8% for the 3-, 4-, and 5-Å channels, respectively.