INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUDING REMARKS REFERENCES

Theoretical studies of biological ion channels have been hampered by a lack of detailed structural knowledge. Until recently, the exact shape of any channel and the positions, densities, and types of dipoles and charge moieties on the protein wall remained unknown. These details are needed to compute the intermolecular potential operating between water molecules, ions, and the protein wall, which is the essential ingredient for theoretical studies of channels using molecular dynamics and, to a lesser extent, Brownian dynamics calculations. A recent report on the crystal structure of the potassium channel (Doyle et al., 1998) has prompted us to carry out electrostatic calculations and simulate the behavior of ions and water molecules in and near the channel to gain an insight into the mechanisms underlying ion permeation and to deduce some of its macroscopically observable properties. There is a need to develop models that can relate the structural parameters of channels to experimental data and thereby build a theoretical framework that can explain different sets of observations. The theoretical description of the potassium channel we give here is produced in the hope of furthering this aim.

The potassium channel is modeled here as a transmembrane lumen, the shape of which corresponds closely to that reported by Doyle et al. (1998), with cylindrical reservoirs containing potassium and chloride ions placed at each end of the channel. Using this basic model, we have examined several key issues from three different perspectivesmacroscopic, semimicroscopic, and microscopic. First, the electrostatic forces experienced by potassium ions at fixed positions in the channel are calculated by using macroscopic approximations. Here, the channel is viewed as a structureless wall made of low-dielectric-strength material, and the water-protein interface is treated as a sharp boundary (Levitt, 1978a, b; Jordan, 1981, 1982, 1983). Second, the trajectories of ions in water interacting with a dielectric boundary are traced using Brownian dynamics simulations. In these simulations, water is treated as a continuum in which ions move under the influence of electrostatic forces and random collisions (Jakobsson and Chiu, 1987; Chiu and Jakobsson, 1989; Bek and Jakobsson, 1994; Li et al., 1998; Chung et al., 1998). Finally, to understand what structural features of the channel render it selectively permeable to potassium ions and how the ion-water geometry undergoes a transformation as an ion moves across the narrow conduit, we have carried out molecular dynamics simulations for all particles in the selectivity filter, as in previous studies on various model pores (see, for example, Roux and Karplus, 1991a; Sansom et al., 1996; Singh et al., 1996; Sankararamakrishnan et al., 1996; Tieleman and Berendsen, 1998). The results of these molecular dynamics calculations are reported in the comparison paper (Allen et al., 1999a).

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUDING REMARKS REFERENCES

The channel model

The transverse section of a model channel, shown in Fig. 1 A, is generated by rotating the two curves (Fig. 1 B) around the symmetry axis (z axis) by 180°. The channel extends from z = 25 Å to 25 Å, with a narrow selectivity filter of radius 1.5 Å and length 12 Å and a wider segment of length 23 Å. The selectivity filter extends toward the extracellular space, whereas the wider pore, whose radius tapers off gradually, extends inward, toward the intracellular space. The radius at the entrance of the channel from the intracellular face is 3 Å. The total interior volume of the channel is 1440 Å³. A cylindrical reservoir of 30 Å radius and variable length is connected to each end of the channel.

View larger version (47K): [in this window] [in a new window]   FIGURE 1   Idealized potassium channel. A transverse section of the model channel shown in A is generated by rotating the curves outlined in B along the symmetry z axis by 180°. The positions of dipoles on the channel wall are indicated in B: , eight of the 16 carbonyl oxygen atoms; , N-terminals of the helix dipoles; , the mouth dipoles.

To investigate how the permeation of ions across the channel is influenced by the presence of fixed charges, we place sets of dipoles in the protein wall with fourfold symmetry around the z axis. First, four rings of four carbonyl groups are placed along the selectivity filter, located at z = 10, 13.33, 16.67, and 20 Å. The negative pole of each carbonyl group (filled circles in Fig. 1 B) is placed 1 Å from the boundary and the positive pole 1.2 Å away from the negative pole, with their orientations perpendicular to the z axis. Second, four helix macrodipoles (open circles), with their N-terminals pointing at the oval chamber near the middle of the channel, are placed 90° apart. The positions of the N-terminals of the helix dipoles are z = 10.66 Å and r = 5.66 Å, and those of the C-terminals are z = 22 Å and r = 17 Å (the length of the dipole is 16 Å). Third, at each entrance of the channel, four "mouth" dipoles (filled diamonds), 5 Å in length, are placed. These are located at z = 22.83 Å and z = 20 Å. In one series of simulations, the strengths of the four helix macrodipoles, 16 carbonyl groups, and eight mouth dipoles are systematically changed to ascertain the strengths that maximize the transfer of ions. In all subsequent series of simulations, the dipole moments of each carbonyl group and each mouth dipole are kept constant at 7.2 and 30 × 10³⁰ Cm, respectively. At each pole of the helix dipoles, we place a charge of ±0.6 × 10¹⁹ C.

Solution of Poisson's equation

For a given configuration of ions and fixed charges in the system, represented by the charge density , the electric potential  is determined by solving Poisson's equation,

(1)

Because the dielectric constant  has different values on either side of the channel boundary (Fig. 1), the solutions of Eq. 1 are subject to the boundary conditions

(2)

where the subscripts 1 and 2 refer to the outside (water) and inside (protein) of the boundary,  is potential, and is the unit vector normal to the surface. Equations 1 and 2 can be solved analytically only for a few channel shapes (Kuyucak et al., 1998), and for the model channel shown in Fig. 1, they have to be solved numerically. This can be achieved with an iterative technique as described by Levitt (1978a) and Hoyles et al. (1996, 1998b). Because the computation of the electric field and potential plays a central role in the ion transport problem, we give a brief description of the technique here and refer to the original references for further details.

The boundary is divided into small sectors of area S_i, each represented by a point charge q_i at its center. The size of each sector varies from 0.3 Å² in places of high curvature to ~14 Ų in flat regions. A total of ~18,000 sectors have been used in the present potassium channel calculations. The boundary conditions (Eq. 2) can be manipulated to relate the surface charge density on each sector to the external electric field E_ex that arises from all of the charges in the system except those in the sector

(3)

where E_ex ·  is determined from the normal derivative of the external potential

(4)

Starting with an initial surface charge density of ⁰(r') = 0, one estimates the potential at the boundary from Eq. 4. This potential is then fed into Eq. 3, and a new density ¹(r') is obtained. Equations 3 and 4 are iterated until the results converge, that is, the difference in  between two successive iterations is less than 0.01% on all sectors. A small error arising from the assumption that the sector is flat is corrected by using the procedure described by Hoyles et al. (1998b). We use a fixed value ₂ = 2 for the protein but vary ₁ in the simulations, which will be denoted simply by .

Lookup tables for electric potentials and fields

The numerical solution of Poisson's equation described above takes a relatively short time on a supercomputer. However, when it has to be repeated at each time step during a computer simulation (typically, many millions of times), the computational cost becomes prohibitively high. We have circumvented this problem by exploiting the huge storage capacity of supercomputers to construct lookup tables for the electric potentials and fields. The method is described and validated by Hoyles et al. (1998a) and used in Brownian dynamics simulations of vestibular channels by Chung et al. (1998). We refer to these references for details of the method and summarize only its essential points here.

The method involves precalculating the electric potential and field on a grid of points for various configurations and storing the results in a number of lookup tables. During simulations, the potential and field at desired points are reconstructed by interpolating between the table entries. Using the superposition principle, it is easy to see that the one- and two-ion configurations are sufficient to construct the potentials in the multiion case. To this end, we break the total electric potential V_i experienced by an ion i into four pieces:

(5)

where the sum over j runs over all of the ions in the system (including the ones in the reservoirs). The four terms in Eq. 5 represent the following components of the total potential:

(i) V_X,i is the external potential due to the applied field, fixed charges in the protein wall, and charges induced by these. Because these quantities do not change during a simulation period, V_X,i depends only on the position of the ion and can be stored in a three-dimensional table.

(ii) V_S,i is the self-potential due to the surface charges induced by the ion i on the channel boundary. Because of the axial symmetry of the channel, V_S,i is independent of the azimuthal angle  and requires only a two-dimensional (2D) table.

(iii) V_I,ij is the image potential due to the charges induced by the ion j. Again because of the axial symmetry, V_I,ij depends only on the relative angle _ij between the two ions and hence can be stored in a five-dimensional (5D) table (instead of six).

(iv) V_C,ij is the Coulomb potential due to the ion j, which is computed directly from

(6)

where r_i and r_j are the positions of the ions.

The electric field is calculated from the derivative of the potential at the grid points and decomposed in exactly the same way:

(7)

with each field component being defined as in the potential (Eq. 5). The only difference is that three values corresponding to the three Cartesian components of the field are stored at each table entry instead of one scalar value in the potential. The dimensions of the 2D, 3D, and 5D tables are 37 × 97, 10 × 171 × 40, and 7 × 119 × 7 × 119 × 14, respectively.

The grid points for the tables are evenly spaced in generalized cylindrical coordinates. This means that, in cylindrical coordinates, the spacing along the z axis (z) is fixed, the angular spacing of points around the z axis () is fixed, and the spacing of points along the radii (r) depends on the radius of the channel. Different tables have different spacings, to minimize the interpolation error for their particular tasks. For example, the 2D table needs very fine radial spacing to minimize error in image repulsion from the channel walls, while the 3D table needs fine linear and angular spacing to accurately represent the field from fixed charges.

The linear spacings (z) are 1.69 Å for the 2D table, 0.96 Å for the 3D table, and 1.37 Å for the 5D table. The angular spacings () are 9.2° for the 3D table and 13.8° for the 5D table. The narrowest part of the channel, the selectivity filter, has a radius of 1.5 Å. The radial spacings (r) here are 0.026 Å for the 2D table, 0.11 Å for the 3D table, and 0.16 Å for the 5D table. The spacings in the widest part of the channel, the cavity of radius 5 Å, are 0.13 Å for the 2D table, 0.51 Å for the 3D table, and 0.76 Å for the 5D table. The spacings in the reservoirs are much larger, as these have a radius of 30 Å.

Born energy

The value of the dielectric constant of water inside the pore is an important open question in channel studies that is often brushed aside by adopting the bulk value of 80. There is no direct experimental information on this quantity, but recent molecular dynamics simulations suggest that it could be much lower than the bulk value in channels with small radii, like the potassium channel (Sansom et al., 1997). Therefore, we take a more flexible approach here and treat  as a variable to be determined from simulations of conductance. The choice of  < 80 in the pore, in turn, raises the question of how to describe the change in  from the bulk value in the reservoir to the lower value in the channel interior. Ideally, one should use a switching function that changes smoothly from one value of  to the other over a given range. However, solution of Poisson's equation with a space-dependent  is a rather complicated computational problem that cannot be tackled with the present numerical techniques. Simplification to a sharp boundary at the channel entrance (similar to the protein boundary) allows solution of the problem by known techniques, but the solutions suffer from instabilities as ions cross this boundary. This problem does not arise with the protein boundary because ions never cross it. As a compromise, we use the same low value of  in the pore and the reservoirs but incorporate the neglected energy difference, which is approximately given by the Born energy,

(8)

as a potential barrier at either entrance of the channel. Here R_B = 1.93 Å, as determined from the enthalpy of hydration for K⁺ ions (Bockris and Reddy, 1970). To avoid complications arising from sharp potential energy changes in Brownian dynamics simulations, we implement this barrier, using a smooth switching function,

(9)

which has continuous first and second derivatives and rises gradually from 0 to E_B as s changes from 1 to 1. Here, z_c = ±22.5 Å is the location of the center of the profile and z = 1.5 Å is its half-width. To give an indication of the barrier heights involved, we note that for  = 20, 40, and 60, E_B = 5.4 kT, 1.8, and 0.6 kT, respectively.

Electrostatic calculations

The potential profile of an ion along the z axis is constructed by solving Poisson's equation with the ion fixed at a given position on the z axis and repeating this procedure at 1-Å intervals. The force experienced by an ion is calculated from the gradient of the potential energy. As will be shown later, the potassium channel is usually occupied by two ions. To visualize the shape of the energy barrier an ion encounters as it attempts to enter a channel that is occupied by one or more ions, we have constructed multiion energy profiles. We move one of the ions from the intracellular space into the channel in 1-Å steps, holding it fixed at each step. We then allow the other ions in the selectivity filter to adjust their positions so that the force on them will be zero, thus minimizing the total energy of the system. The minimization is performed at each step, and the positions of the ions and the total energy are recorded. This corresponds to the total electrostatic energy required to bring in the charge on the ions from an infinite distance in infinitesimal amounts, and it is given by

(10)

where the indices i and j range over all of the ions. The potential terms in Eq. 10 assume the same significance as in Eqs. 5 and 9. The factors of 1/2 in the middle three terms arise from the integration of charge during build-up. Using a modified version of the steepest descent algorithm (see Press et al., 1989), the total energy of a multiply occupied channel given in Eq. 10 is successively minimized until the forces on the free ions vanish.

Brownian dynamics

Brownian dynamics (BD) simulations are used to predict the channel conductance under various conditions and to deduce the optimum strength of various dipole groups for the most efficient transfer of ions across the membrane. In BD, the motion of an ion with mass m_i and charge q_i in a fluid is governed by the Langevin equation,

(11)

By solving the Langevin equation at discrete time steps for each ion in the system, for each Cartesian component (x, y, z) of the velocity, we determine the number of ions crossing the channel in a fixed simulation period. The three terms on the right-hand side of Eq. 11 correspond to, respectively, an average frictional force with a friction coefficient given by m_ii, the stochastic force arising from random collisions with water molecules, and the total systematic force acting on the particle. The systematic force is composed of short- and long-range forces. As outlined above, the latter are obtained from the numerical solution of Poisson's equation on a grid of points and stored in a number of lookup tables. During simulations, the field and potential at desired positions are reconstructed by interpolating between the table entries. The short-range forces include the Born energy barrier introduced in Eqs. 8 and 9, and part of the ion-ion interaction U_II and the ion-wall interaction U_IW. During the simulations, ions come closer to each other at times than the sum of their radii. This activates a short-range repulsion arising from the overlap of their electron clouds, given by

(12)

where r₁₂ is the ion-ion distance; R_i, i = 1, 2, are the Pauling radii of ions; and F₀ is the magnitude of the short-range force at contact. For the ion-protein wall potential U_IW, a similar form is used:

(13)

where R_W is the radius of the atoms making up the wall and (z) = R_C(z)  r is the distance of the ion from the channel wall at R_C(z). We use R_W = 1.4 Å and F₀ = 2 × 10¹⁰ N in both Eqs. 12 and 13, which is estimated from the ST2 water model used in molecular dynamics (Stillinger and Rahman, 1974).

To simulate the effects of short-range forces more accurately, we use a multiple time-step algorithm in our BD code. A shorter time step of 2 fs is used in the mouth regions of the channel, where U_B is active, and in the narrow regions, where U_IW is expected to contribute significantly. A long time step of 100 fs is used elsewhere. Specifically, there are two short time step bands, 25 < z < 15 and 7.5 < z < 25, comprising both entrances and the selectivity filter. If an ion is in one of these bands at the beginning of a 100-fs period, it is simulated by 50 short steps instead of one long step; so synchronization between the ions is maintained. Long-range forces are calculated normally at the start of the 100-fs period and are assumed to be constant throughout. The ion-ion interactions are normally treated using the long time steps, except when both ions are in one of the above bands.

Technical details of simulations

We solve the Langevin equation using the BD algorithm devised by van Gunsteren et al. (1981, 1982), which consists of the following computational steps:

Step 1. Compute the electric force F(t) = q_iE_i acting on each ion at time t_n and calculate its derivative [F(t_n)  F(t_n1)]/t.

Step 2. Compute a net stochastic force impinging on each ion over a time period of t from a sampled value of F_R(t).

Step 3. Determine the position of each ion at time t_n + t and its velocity at time t_n by substituting F(t_n), its derivative F'(t_n), and F_R(t) into the solutions of the Langevin equation (Eqs. A6 and A7 of Hoyles et al., 1998a).

Step 4. Repeat the above steps for the desired simulation period.

Simulations under various conditions, each lasting for 1,000,000 time steps (0.1 µs), are repeated many times, mostly for 5 to 10 trials. For the first trial, the positions of ions in the reservoirs are assigned randomly with the proviso that the minimum ion-ion distance should be at least 2.7 Å. For successive trials, the positions of the ions in the last time step are used as the initial starting positions of the following trial. The current (given in pA) is determined from the total number of ions traversing the channel over the simulation period.

Fixed numbers of potassium and chloride ions are placed in each reservoir, and the height of the cylindrical reservoir is adjusted to give a desired ionic concentration. As ions are forbidden to approach the wall of the reservoir within 1 Å, the effective radius of the cylindrical reservoir is 29 Å. For example, if 13 sodium and 13 chloride ions are placed in each reservoir and the desired ionic concentration is 300 mM, the height of the cylindrical reservoirs is adjusted to 27 Å.

When the ionic concentration in the reservoirs is high, ions at times are able to jump large distances and end up very close to another ion. The forces at the next time step in such instances can be very large, and the affected ions may leave the system. To correct this problem, we check ion-ion distances at each time step. If two ions are within a "safe distance" (see Chung et al., 1998), chosen to be 3/4 of the sum of the ionic radii (Pauling, 1942), then their trajectories are traced backward in time until such a distance is exceeded. By performing these checks and corrections, the system is well behaved over the simulation, even for very high concentrations. Such a minor readjustment of the position of an ion is needed about once every 100 time steps when the reservoirs and the channel contain 52 ions. The steep repulsive force at the dielectric boundary due to the image charges and the ion-protein potential U_IW is usually sufficient to prevent ions from entering the channel protein. We ensure that no ions appear inside the channel protein by erecting an impermeable hard wall 1 Å from the water-protein interface. Any ion colliding with this wall is elastically scattered. A similar hard wall is implemented for the reservoir boundaries.

To ensure that the desired intracellular and extracellular ion concentrations are maintained throughout the simulation, a stochastic boundary is applied. When an ion crosses the transmembrane segment, an ion of the same species is transplanted so as to maintain the original concentrations on both sides of the membrane. For example, if a potassium ion from the left-hand side of the channel crosses the pore and reaches the imaginary plane at z = 25 Å, then a potassium ion located at the furthermost right-hand side reservoir is taken out and placed in the far left-hand side of the left reservoir. When transplanting ions, we choose a point no closer to another ion than the defined safe distance. The stochastic boundary trigger points, located at z = ±25 Å, are checked at each time step of the simulation.

We represent the potential difference across the channel by an applied electric field of constant strength E. In the absence of any dielectric boundary, the potential difference across a channel of length d would be E × d. The presence of a dielectric boundary and dipoles on the protein wall, however, severely distorts the field. Thus, the precise potential difference across the channel will depend on the selected reference points at the two sides of the potassium channel. For simplicity, we apply a field strength of 10⁷ V m¹ and refer to it as an applied potential of 100 mV. To construct the current-voltage relationships and accurately determine the reversal potentials with different ionic concentrations in the reservoirs, however, we apply a fixed field strength in the presence of all of the dipoles (but without placing any ions in the reservoir) and then measure the electric potential at the middle of each reservoir (usually z = ±40 Å) at the central axis. The current is then plotted against the potential difference between these two reference points.

The BD program is written in Fortran and 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 52 ions in the reservoirs, the CPU time needed to complete a simulation period of 1.0 µs (10 million time steps) is ~28 h.

Throughout we quote energy in room temperature units (kT) and dipole moment in Coulomb-meter (Cm). We note 1 kT = 4.11 × 10²¹ J or 0.592 kcal/mol and 1 Debye = 3.33 × 10³⁰ Cm. The following physical constants were employed in our calculations (note that the friction coefficient is related to the diffusion coefficient via the Einstein relation D = kT/m):

Masses: m_K = 6.5 × 10²⁶ kg, m_Cl = 5.9 × 10²⁶ kg.

Diffusion coefficients: D_K = 1.96 × 10⁹ m²/s, D_Cl = 2.03 × 10⁹ m²/s.

Ionic radii: R_K = 1.33 Å, R_Cl = 1.81 Å.

Room temperature: T = 298 K.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION CONCLUDING REMARKS REFERENCES

Dipoles and energy profiles

As an ion approaches the boundary between an aqueous solution and the protein wall, it experiences an electrostatic repulsion due to induced charges at the boundary. In computing the potential energy of an ion as it moves along the central axis, we assume initially that the dielectric constant  in the reservoirs and the pore is 60. The energy of transition from bulk water, estimated from the Born energy, is incorporated as a potential barrier at the channel entrance, as explained in the Methods section.

In the absence of any charge moieties on the protein wall, an ion attempting to traverse the channel encounters a significant energy barrier. The potential energy at a fixed position of an ion is computed numerically, and the calculation is repeated at 1-Å intervals. The profile presented to the ion as it moves from inside (left) to outside (right) increases slowly at first and then rises steeply in the narrow selectivity filter, reaching a peak of 20 kT, as shown in Fig. 2 A (curve labeled a). Four rings of dipoles, with four carbonyl groups in each ring, placed along the selectivity filter, transform a section of the barrier into a well (b), as do four helix dipoles placed just below the selectivity filter (c). As will be shown later, two sets of additional mouth dipoles are needed to render the channel permeable to ions. When all three sets of dipoles16 carbonyl groups, four helix macrodipoles, and eight mouth dipolesare placed along the channel wall, the profile an ion encounters while traversing the central axis of the channel is a deep potential well (d).

View larger version (15K): [in this window] [in a new window]   FIGURE 2   Electrostatic energy profile of a potassium ion traversing the channel. (A) The potential energy of an ion along the central axis is calculated numerically at 1-Å intervals. The profile a is the potential energy seen by the ion, in the absence of any charge moieties on the channel wall; profile b is obtained with four sets of four carbonyl groups, ( in the inset); profile c is obtained with four helix macrodipoles, () in the inset); and profile d is obtained with the carbonyl groups, helix dipoles, and two rings of mouth dipoles, (, in the inset). (B) Potential profiles seen by individual ions while the other is fixed at the equilibrium position (indicated by arrows). The applied field corresponds to a potential difference of ~150 mV between the channel entrances. The first ion is placed at the equilibrium position z = 10 Å, and the profile seen by the second ion is computed (solid line). The second ion is now placed at z = 19.75 Å and the profile seen by the first ion is computed (broken line). (C) The potential profile encountered by a third ion, as it moves from left to right. At each fixed position of the third ion, the stable configuration of the first two ions is determined iteratively, and then the total energy of the assembly is computed. This process is repeated at 1-Å intervals from z = 40 to 10 Å to produce the profile shown as a solid line. Unless stated otherwise, the strengths of each carbonyl group, mouth dipole, and helix dipole are, respectively, 7.2, 30, and 96.3 × 1030 Cm.

The potential well created by the dipoles, reaching a depth of nearly 30 kT, attracts cations. An ion, upon entering the channel, will proceed toward the bottom of this well. A second ion entering the channel sees a different profile, altered by the presence of the first ion. The well in Fig. 2 A (d) is deep enough to hold two ions in a stable configuration. Through an iterative energy minimization procedure, one can determine the equilibrium positions of the pair of ions in the well. The potential profile seen by either ion while the other is fixed at the equilibrium position, in the presence of an applied field of 1.5 × 10⁷ V/m, is shown in Fig. 2 B. At these positions (indicated by arrows), the z-component of the force experienced by the ions is zero. The two-ion potential profiles exhibit relatively deep wells that may attract a third ion. In Fig. 2 C we show the potential profile seen by a third ion moving into the channel from the left. Here the potential is calculated at a given position of the third ion after the equilibrium positions of the first two ions are found. There is a shallow well near the entrance of the channel, produced by the ring of mouth dipoles. Once in the well, the third ion will be delayed until random Brownian motion allows it to escape. We note here that the potential minimum is along the central channel axis, so that ions are preferentially funneled along it. The repulsive force from the induced surface charges swings into action whenever an ion strays from the central axis, pushing it back to the axis. The corresponding electric potential profile along the radial axis is similar in appearance to a harmonic well, except that it rises much more sharply near the boundary (see Hoyles et al., 1996).

The relative contributions of various charge groups in establishing a potential gradient in the channel are summarized in Fig. 3. The curves reveal the electric potential; they differ from the potential energy curves in Fig. 2 A in that there are no induced surface charges due to ions. Four mouth dipoles at each end of the channel produce a potential well (Fig. 3, (a)). The depths of the wells near the intracellular and extracellular entrances reach, respectively, 226 mV and 202 mV. The carbonyl groups lining the selectivity filter produce a steep potential well that reaches a maximum depth of 681 mV at z = 15.8 Å (b). A broader well encompassing nearly the entire extent of the channel is produced by the helix dipoles (c). It reaches a minimum of 564 mV at z = 10 Å. The potentials produced by these three groups of dipoles sum algebraically when all of the dipoles are placed, as shown in (d). At z = 14 Å, the potential experienced by the test charge reaches a minimum of 1250 mV.

View larger version (14K): [in this window] [in a new window]   FIGURE 3   Electric potential generated by dipoles. The electric potential is calculated numerically at 1-Å intervals. The profiles labeled a, b, and c are obtained, respectively, with eight mouth dipoles, 16 carbonyl groups, and four helix dipoles. The potential profile, labeled d, resulting from all three groups of dipoles, is the sum of a, b, and c.

From these electrostatic calculations, we deduce that the channel is normally occupied by two cations. Conduction is unlikely to take place unless these ions are resident in the pore to reduce the energy well created by the charge moieties. Moreover, for the channel to conduct ions, the effective dielectric constant needs to be quite large.

Dependence of conductance on dipole strengths

For the channel to conduct ions, there is a narrow range of moments various dipole groups must possess. Using BD simulations, we have determined how the magnitude of currents flowing across the channel varies with dipole strengths and the effective dielectric constant in the channel lumen.

As shown in Fig. 4 A, the conductance increases rapidly as the moment of each carbonyl group is increased from 0 to 7.2 × 10³⁰ Cm. The current begins to decline when the moment is further increased to 14.4 × 10³⁰ Cm. The three curves illustrated in Fig. 4 A are obtained by letting the effective dielectric constant of the pore be 80 (top curve), 70 (middle curve), and 50 (bottom curve). The charge placed on the terminals of each helix dipole and the strength of each mouth dipole are kept constant at 0.6 × 10¹⁹ C and 30 × 10³⁰ Cm, respectively. In this and subsequent figures, unless stated otherwise, each point is the average of five simulations, with each simulation period lasting for 100 ns. The error bar accompanying a data point is one standard error of means and is not shown if it is smaller than the size of the data point. Again, unless noted otherwise, 13 potassium and 13 chloride ions are placed in the left-hand reservoir (representing the intracellular space), whose volume is adjusted so as to give an ionic concentration of 300 mM, and the same number of ions is placed in the right-hand reservoir (representing the extracellular space). The applied electric field between the two ends of the reservoirs produces a potential difference of ~150 mV, inside positive with respect to outside. The peak current is always obtained when the strength of the carbonyl groups is fixed at 7.2 × 10³⁰ Cm (2.16 Debye), irrespective of the assumed dielectric constant. In Fig. 4 B, the variation of currents with the dipole moment is determined at three different applied potentials, 150 mV, 200 mV, and 250 mV, while keeping  = 60 throughout. Again the current peaks at about the same dipole strength.

View larger version (19K): [in this window] [in a new window]   FIGURE 4   Changes in channel conductance with the strength of carbonyl groups. Simulations under various conditions, each lasting 100 ns, are repeated five times. The height of each reservoir is adjusted to give an ionic concentration of 300 mM, and inside is made 150 mV positive with respect to the outside. The current (in pA) is determined from the total number of ions traversing the channel over the simulation period of 0.5 µs. (A) Keeping the moments of each helix dipole and mouth dipole constant at 96.3 and 30 × 1030 Cm, the strength of each carbonyl group is systematically changed. The current (in pA) is plotted against dipole strength for three values of the dielectric constant,  = 50, 70, and 80. (B) The current is plotted against the strength of carbonyl groups for three values of applied potential, 150 mV, 200 mV, and 250 mV, while keeping  constant at 60.

The results of simulations showing the variation of current with the strengths of mouth dipoles (A) and helix dipoles (B) are illustrated in Fig. 5. The dipole moment of each carbonyl group in the selectivity filter is kept at 7.2 × 10³⁰ Cm throughout. The current flowing across the channel is largest when the charge on each of the four helix dipoles is 0.6 × 10¹⁹ C. Similarly, the current is maximum when the strength of each of the mouth dipoles is 30 × 10³⁰ Cm. Fig. 5 reveals, as does the previous figure, that the dielectric constant of the channel has a pronounced effect on the permeability of the channel. With optimum pore helix and mouth dipole strengths,  = 60 gives a physiological conductance of ~40 pS, as found experimentally (Schrempf et al., 1995). The channel conductance is progressively suppressed when the dielectric constant in the pore is lowered and no conduction takes place, with a driving force of 150 mV, when   40. 

View larger version (19K): [in this window] [in a new window]   FIGURE 5   Changes in channel conductance with the strength of mouth dipoles (A) and helix dipoles (B). (A) The strength of each mouth dipole is changed systematically while keeping the moment of each carbonyl group and helix dipole constant at 7.2 and 96.3 × 1030 Cm. The current is plotted against the strength of mouth dipoles for three different values of dielectric constant,  = 60, 70, and 80. (B) The charge placed on each helix dipole is changed systematically while keeping the moments of each carbonyl group and mouth dipole constant at 7.2 and 30 × 1030 Cm. The current (in pA) is plotted against dipole strength for three values of the dielectric constant,  = 50, 60, and 70.

Here and in all subsequent series of simulations, we assume that the channel possesses the strengths of various dipole groups, which enable the maximum number of ions to be translocated across the channel for a given driving force, that is, the dipole strengths for each of eight mouth dipoles, 16 carbonyl groups, and four helix dipoles are, respectively, 30, 7.2, and 96.3 × 10³⁰ Cm.

Effects of dielectric constant and diffusion coefficient on currents

From the results given in the previous section, it is clear that, for the channel to conduct, the effective dielectric constant  in the pore must be high. In other words, water molecules resident in the pore must not be tightly bound to the protein but be able to rotate relatively freely so as to reduce the interaction energy between the ions and the charges located on the channel wall. In Fig. 6, we examine further the influence of  on the magnitude of current flowing across the channel. The depth of the energy well created by dipoles increases as  is lowered. An example of the energy well created by four mouth dipoles located near the channel entrance, when there are two ions resident in the selectivity filter (c.f., Fig. 2 C), is illustrated in Fig. 6 A. Here,  is assumed to be 30, and a potential of 300 mV is applied across the channel. An ion attempting to cross this well encounters a barrier V_B, the height of which decreases monotonically with increasing , as shown in Fig. 6 B. Increasing the applied potential from 150 mV to 300 mV reduces the barrier height by ~1.5 kT. A steep increase in the barrier height as  is lowered suggests that the channel will not conduct ions if  in the pore is less than 40. The inference drawn from electrostatic calculations is in accord with the results obtained from BD simulations. The current across the channel under the driving force of 150 mV, 200 mV, 250 mV, and 300 mV is plotted against  in Fig. 6 C. These four curves broadly mirror the way the barrier height increases with . The current ceases to flow when the barrier height reaches 7 kT.

View larger version (16K): [in this window] [in a new window]   FIGURE 6   Effects of the effective dielectric constant on conductance. (A) An energy well near the channel entrance created by four mouth dipoles is encountered by an ion, given that there are two resident ions in the selectivity filter. The energy minimization is carried out with  = 30 and an applied potential of 300 mV. An ion upon entering the well must surmount a barrier of height VB to traverse the channel. (B) The height of the barrier VB is plotted against the effective dielectric constant for two different values of applied potential, 150 mV () and 300 mV (). (C) The inward currents under the driving force of 300 mV (), 250 mV (), 200 mV (), and 150 mV () are plotted against the effective dielectric constant. Each point is derived from a simulation period of either 1 µs (150, 200, and 250 mV) or 2 µs (300 mV).

The diffusion coefficient of potassium ions D_K in bulk electrolyte solutions is 1.96 × 10⁹ m²/s. This value is reduced when an ion is diffusing through a narrow tube (Roux and Karplus, 1991b; Chiu et al., 1993; Lynden-Bell and Rasaiah, 1996; Smith and Sansom, 1997; Allen et al., 1999b). The magnitude of the diffusion coefficient of an ionic species depends on, among other things, the radius of the cylinder and the composition of the wall. In the accompanying paper (Allen et al., 1999a), we show that D_K in the wider segment of the potassium channel, including the oval chamber, is nearly the same as that in bulk solutions, whereas that in the selectivity filter is on average of the bulk value. The following series of simulations are carried out to assess how much the channel conductance is suppressed by a low D_K in the narrow filter.

When ions enter the channel segment extending from z = 7.5 to z = 25 Å, their motions are determined by a D_K that is lower than the bulk value. Fig. 7 shows the current across the channel as a function of D_K at three different values of dielectric constants,  = 60 (A) and  = 50 and 70 (B). The filled and open circles in Fig. 7 A represent, respectively, the outward and inward currents. The applied potential across the channel and the ionic concentration in the reservoir are kept constant at 200 mV and 300 mM, respectively. In contrast to bulk conductance, where current is proportional to the diffusion coefficient, the current in the potassium channel depends on D_K in a nonlinear fashion. It decreases with decreasing D_K at a very slow rate at first (until D_K is reduced to ~0.5 of its bulk value) and then becomes more or less proportional to D_K. When D_K is reduced to of the bulk value, the current is only suppressed by ~30%.

View larger version (14K): [in this window] [in a new window]   FIGURE 7   Effects of the diffusion coefficient on conductance. The diffusion coefficient of potassium ions in the selectivity filter, extending from z = 7.5 to 25 Å, is progressively lowered from the bulk value of 1.96 × 109 m2/s, and its effect on conductance is examined. (A) The outward () and inward () current is plotted against diffusion coefficient in the selectivity filter, with  = 60 and an applied potential of 200 mV. (B) The outward current is plotted against diffusion coefficient in the selectivity filter for  = 70 (upper curve) and 50 (lower curve). Each point in A and B is derived from a simulation period of 1 µs.

Current-voltage relationships

The current-voltage relationships, shown in Fig. 8, are obtained with symmetrical solutions of 300 mM in both reservoirs. The diffusion coefficient in the selectivity filter is assumed to be the same as that in bulk electrolytes. Because the effective dielectric constant  of the channel is unknown, we have determined the current-voltage curves, assuming  to be 60, 70, and 80. The curves derived from these three conditions all reveal several distinct features. First, at any given applied potential, the outward current is larger than the inward current. Second, the magnitude of current across the channel at any given driving force increases steadily with increasing dielectric constant. The outward current at 100 mV is 6.7 ± 1.2, 11.8 ± 2.1, and 15.0 ± 1.0 pA when  is assumed to be 60, 70, and 80, respectively (Fig. 8). Because the current begins to saturate with increasing ionic concentrations (see later), the conductance at 150 mM K⁺ will be slightly larger than 33, 59, and 75 pS at these three values of dielectric constants. Third, the relationship is approximately linear when the applied potential is less than 100 mV, but it deviates 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 an ion when the driving force is large. Thus, in the presence of a barrier, the ohmic current-voltage relationship will be modified by a function of the form

(14)

where  is the limiting conductance at large V,  is a dimensionless parameter, and V_B1 and V_B2 are the right and left barrier heights. The justification for fitting the data with this equation is given by Chung et al. (1998). Rectification arises from the fact that ions moving into and out of the cell see different barriers. The solid lines fitted through the data points (Fig. 8, A, B, and C) are calculated from Eq. 14.

View larger version (10K): [in this window] [in a new window]   FIGURE 8   The current-voltage relationships. The current measured at various applied potentials is obtained with symmetrical solutions of 300 mM in both reservoirs. The solid lines fitted through data points are calculated from Eq. 14. The values of  used for A, B, and C are 60, 70, and 80, respectively.

The current-voltage relationships obtained with asymmetrical ionic solutions in the two reservoirs are shown in Fig. 9. The curves exhibited in the figure are obtained by assuming that  in the channel is, respectively, 60 (A) and 70 (B). The ionic concentrations inside and outside are 500 mM and 100 mM. The solid lines fitted through the data points are obtained from Eq. 14, multiplied by the Goldman factor of the form

(15)

As expected, the asymmetry between the inward and outward currents is accentuated. As with the symmetrical solutions, an increase of  from 60 to 70 causes an increase in the magnitude of currents flowing across the channel. The zero current of the two current-voltage relationships appears to be somewhere between 25 mV and 50 mV.

View larger version (11K): [in this window] [in a new window]   FIGURE 9   The current-voltage relationships. The current measured at various applied potentials is obtained with asymmetrical solutions for  = 60 (A) and 70 (B). The concentration in the reservoir representing the intracellular space is 500 mM, whereas that representing the extracellular space is 100 mM. The solid lines drawn through the data are calculated with Eq. 14 multiplied by the Goldman factor. For C, obtained with  = 80, the measured concentrations in the intracellular and extracellular reservoirs are 482.0 mM and 71.5 mM () and 385.3 mM and 176.2 mM (). The open downward arrows indicate the reversal potential calculated from the Nernst equation. The predicted Nernst potentials taking the activity coefficients into account are indicated with filled downward arrows. The values of activity coefficients used for 71.5, 176.2, 386.3, and 482.0 mM are, respectively, 0.79, 0.73, 0.67, and 0.63.

To ascertain how closely the measured reversal potentials match those predicted by the Nernst equation, we estimate currents flowing across the channel with two different ionic concentrations in the reservoirs and under various applied potentials. The concentrations of K⁺ in the extracellular and intracellular aspects of the channel are computed from the average number of ions in the reservoirs throughout the simulation periods. The measured ionic concentrations in the left and right reservoirs in one series of simulations are 71.5 and 482.0 mM, and in another series of simulations are 176.2 and 385.3 mM. Fig. 9 C shows the currents flowing across the channel at various applied potentials. Because the net current for these driving forces is small, the total simulation period of 3 µs is used to derive each data point. For the same reason, we use  = 80 for the effective dielectric constant of the channel, which results in a larger current flow. The reversal potential for each asymmetrical solution is estimated by fitting a polynomial curve through the data points (solid line in Fig. 9 C). There are small but consistent discrepancies between the reversal potentials deduced from simulations and those predicted from the Nernst equation (indicated with open downward arrows). The zero currents occur at 45 mV and 17 mV when the concentration ratios in the two reservoirs are, respectively, 6.7:1 and 2.2:1. The predicted reversal potentials are 48.1 mV and 19.7 mV. These discrepancies between the predicted and measured zero currents disappear if we take the activity coefficients of KCl at the measured ionic concentrations into account (Weast, 1983), as indicated by the filled arrows in Fig. 9 C. From a number of I-V curves obtained with asymmetrical solutions, we conclude that the zero current occurs at a potential predicted by the Nernst equation within the errors of simulations.

Ions in the channel

It is of interest to note where in the channel ions dwell predominantly. To compute the average number of ions inside the channel, we divide the channel into 32 thin sections and compute the time averages of potassium ions in each section. When a potential of 200 mV is applied so as to produce an outward current, two ions on average tend to reside in the channel. The preferred positions where ions dwell are in the selectivity filter at z = 9.4, 14.1, and 23.4 Å, as shown in Fig. 10 A. We note here that, although the histogram (Fig. 10 A) shows three distinct peaks near the selectivity filter, there are on average 1.5 ions in this region, as can be deduced by summing the heights of the bars. A similar sum for the peak near the intracellular entrance gives 0.5 ions; that is, an ion is present there 50% of the time. The preferential positions of the ions in the channel are shifted when the direction of the current is reversed by making the inside negative with respect to the outside. Under this condition, two ions mainly linger around z = 9.4 and 17.2 Å (Fig. 10 B). Thus the preferred locations of ions in the channel depend on, among other factors, the direction and the strength of the applied field.

View larger version (19K): [in this window] [in a new window]   FIGURE 10   Concentrations of potassium ions in the channel. The channel is divided into 32 sections, as indicated in the inset, and the probability of ions present in each section over a simulation period of 0.5 µs is tabulated (bars). The ionic concentration in the reservoirs is 300 mM. The applied electric field in A is 2 × 107 V/m, such that the inside (left-hand side) is ~200 mV positive with respect to outside (right-hand side). The direction of the field is reversed to obtain the distribution shown in B.

To better illustrate the behavior of ions under the influence of the electric and stochastic forces, we bisect the channel and denote the number of ions on the left-hand and right-hand sides by [n_l, n_r]. The occupation probabilities of distinct states are tabulated in Table 1 for five different potentials. At the bottom of the table, we give the average number of ions resident in the channel, which is about two, regardless of the applied potential. In view of the rapid change in occupation probabilities of different states (Table 1), this appears as quite remarkable, and reinforces the earlier inferences made from electrostatic calculations that two ions must be resident in the pore for conduction to take place.

View this table: [in this window] [in a new window]   TABLE 1   Occupation probabilities of multiion states [nl, nr] and the average number of ions in the channel nl + nr for different applied potentials

When the applied potential is 100 mV, which is relevant for the operation of the potassium channel, the most common state is the [0, 2] state. That is, no ion is present in the first half (intracellular side) of the channel, while two ions are present in the second half (extracellular side), as illustrated schematically in the upper panel of Fig. 11 A. In addition to the states listed in Table 1, there are five other distinct states that are observed during the total simulation period of 0.5 µs (or 5 million time steps), but the frequencies of their occurrences are less than 1%. About 32,000 transitions occur between these states when the snapshot of the channel state is taken once every picosecond. The most common transitions are between [0, 2] and [0, 1], and between [1, 2] and [1, 1], which corresponds to the process: driven by thermal energies, one of the two ions in the second half of the channel escapes and then reenters. The forward and backward transitions between these two sets of states account for 64% of the total transitions. Less frequent transitions (~20% of all transitions) are between [0, 2] and [1, 2]. Finally, transitions between [0, 1] and [1, 1] account for 6% of the total transitions, while the forward transition between [1, 1] to [0, 2] accounts for only 0.3% of the total transitions.

View larger version (13K): [in this window] [in a new window]   FIGURE 11   Ions in the channel. (A) With an applied potential of 100 mV, the most commonly occurring configuration is exhibited in the upper panel. Two ions in the selectivity filter are indicated as open circles. The mouth dipoles create an energy well, whose depth VB is ~4.4 kT (middle panel). The average time it takes for an ion to enter the well and then to climb out of the well (t1) and the time to reach the selectivity filter (t2) are indicated in the lower panel. (B) The average time it takes for an ion outside the channel to reach the deepest section of the energy well (z = 20 Å) after the previous ion successfully climbed out of the well () and the mean velocity of an ion between z = 10 and z = +14 Å () are plotted against the applied potential. The concentrations of the reservoirs are kept constant at 300 mM.

For an ion to traverse the channel from inside to outside, one of the rare sequences of state transitions that must take place is [0, 2]  [1, 2]  [0, 3]  [0, 2]. With an applied potential of 100 mV, ~40 ions traverse the channel in 1 µs, corresponding to a processing time of 25 ns per ion. Analysis of the trajectories of ions in the assembly reveals that the rate-limiting step in conduction is the time it takes for a third ion to stumble into the channel entrance and then climb up the energy barrier displayed in the middle panel of Fig. 11 A. This time, indicated by t₁ in the lower panel of Fig. 11 A, is almost equal to the total processing time. The drift velocity of the ion that successfully climbs over the barrier is about one order of magnitude larger than the reservoir value. Thus it takes only t₂ = 1.0 ns for the ion to go from the barrier to the selectivity filter (z = 10 to 14 Å), which is much shorter than one would naively expect from diffusive kinematics. The ejection of the rightmost ion in the selectivity filter is concurrent with the process t₂; hence for all practical purposes, t₁ determines the processing time.

Further analysis of the time t₁ in terms of the channel access and barrier transit times is complicated because the transition region between these two processes is not well defined. In the following, we use the most apparent choice, namely, the imaginary plane at the channel entrance (z = 25 Å), to study the voltage and concentration dependence of the individual waiting time. In Fig. 11 B, we plot, as a function of the applied membrane potential, the average time it takes for an ion to reach the bottom of the well (z = 20 Å) after the previous ion has successfully climbed out of the well. This waiting time, or access time, decreases steadily from 18.5 ns at the applied potential of 50 mV to 1.45 ns at 300 mV (open circles). Also shown is the average speed of the ion (to travel from z = 10 to +14 Å) after it has successfully climbed over the barrier (filled circles). The drift velocity in the channel, as expected, increases mono