INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

In the previous article (Moy et al., 2000, hereafter called article I), we have tested the validity of the mean-field approximation in Poisson-Boltzmann (PB) theory, which is commonly used in potential energy calculations in ion channels. The PB theory is limited to equilibrium situations; to describe non-equilibrium processes such as ion transport, another continuum theory that is widely known as the Nernst-Planck (NP) electrodiffusion equation is used. The NP equation combines Ohm's law for drift of ions in a potential gradient with Fick's law of diffusion due to a concentration gradient (hence the name "drift-diffusion equation" is used in some fields). When the potential in the NP equation is determined from Poisson's equation in a self-consistent manner, the combined system of equations form the Poisson-Nernst-Planck (PNP) theory, which provides a premium description of ion transport problems in many branches of physics and chemistry (e.g., Ashcroft and Mermin, 1976; Bockris and Reddy, 1970; Mason and McDaniel, 1988; Newman, 1991; Weiss, 1996). As in the case of the PB theory, these applications usually involve bulk conditions with system sizes much larger than the Debye length, and the validity of the underlying mean-field approximation is well established. Recent applications of the NP and PNP theories in ion channels (see Levitt, 1986; Cooper et al., 1988; Hille, 1992; Eisenberg, 1996, 1999 for reviews and further references), in contrast, involve systems with rather few ions and with dimensions smaller than the Debye length. Under these conditions, one would intuitively expect that keeping the integrity of ions would be essential to gain a realistic physical description of the system, and the validity of the continuum approaches, where ions are represented as a continuous charge density, would be largely compromised. The most direct way of checking the validity of the PNP theory is to compare its predictions for various physical quantities (e.g., current and concentration) with those obtained from Brownian dynamics (BD) simulations, where individual ions are treated explicitly. The importance of such a test of PNP theory has been stressed in a recent series of commentaries on ion permeation by all participants (Levitt, 1999; McClesky, 1999; Miller, 1999; Nonner et al., 1999). Although molecular dynamics simulations (Roux and Karplus, 1994) are not yet at a stage to replace PNP theory in case of failure, BD simulations currently provide a genuine alternative for studying ion permeation in channels (Li et al., 1998; Chung et al., 1998, 1999; Hoyles et al., 1998a).

In this article we test the PNP theory by comparing its predictions for conductance and concentration profiles in cylindrical channels and a potassium channel with those of BD simulations. We emphasize that both theories are applied to three-dimensional (3-D) channels without any simplifying assumptions that would reduce them to equivalent 1-D problems. Extension of both theories from effective 1-D channels to realistic 3-D cases has been achieved very recently (see Kurnikova et al., 1999 for PNP, and the BD references quoted above). The 3-D aspect of the channel structure is very important in settling such questions as the amount of shielding of dielectric forces on ions. In this respect, the earlier 1-D BD simulations of ion channels (Cooper et al., 1985; Jacobsson and Chiu, 1987; Bek and Jacobsson, 1994) provide only a limited testing ground for the continuum theories.

We note that the continuum description of water in both BD and PNP is strictly valid in bulk situations, and the channel-water interactions are expected to play a role in ion permeation. These interactions can be directly taken into account in molecular dynamics studies. However, the infeasibility of molecular dynamics simulation of ion permeation with the currently available supercomputers necessitates a more phenomenological approach to the problem. One hopes that the effective parameters used in the phenomenological approaches (such as diffusion coefficient and dielectric constant) will all be determined from molecular dynamics studies eventually. This will provide a bridge between the microscopic and macroscopic approaches and a justification for the use of the latter theories. In the mean time, it is important to test the validity of various approximations going into the phenomenological theories to assess their suitability as models of membrane channels. Our comparison of PNP and BD in this article is carried out in this spirit. A summary of some of the results illustrated here was communicated elsewhere (Corry et al., 1999).

	THEORETICAL METHODS

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

PNP theory

The PNP approach to ion permeation in membrane channels has been used in numerous papers in the last decade. Here we give an outline of the theory, and refer to the recent review articles by Eisenberg (1996, 1999) for further details and references. More recent references can be found in Chen et al. (1997, 1999), Nonner and Eisenberg (1998), and Kurnikova et al. (1999).

In continuum theories, the flux J of each ion species is described by the NP equation, which combines the diffusion due to a concentration gradient with that from a potential gradient

(1)

here D, ze, and n are, respectively, the diffusion coefficient, charge, and number density of the ions of species . Note that n (in SI units) is related to the concentration of ions c (in moles/liter) through n = 10³N_Ac. The potential  in Eq. 1 is determined from the solution of Poisson's equation

(2)

where  is the dielectric constant, the sum over  gives the charge density associated with the mobile ions in the electrolyte, and _ex represents all the other external charge sources such as fixed or induced charges on a boundary. In the PNP theory, Eqs. 1 and 2 are solved simultaneously, yielding the potential, concentration, and flux of ions in the system. Note that both the ion concentration and flux are described by continuous quantities corresponding to macroscopic, space-time averages of microscopic motion of individual ions.

Due to their nonlinear nature, the PNP equations are notoriously difficult to solve analytically except in some very special cases, e.g., the classic Goldman-Hodgkin-Katz equation (Hille, 1992). More recent discussions of the analytical treatment of the PNP equations can be found in Syganow and von Kitzing (1995, 1999a, b). Here we consider the basic formalism of the PNP together with some special cases to indicate where and why the PNP theory may break down. These solutions will also be used in checking the accuracy of the numerical results.

When J = 0 in Eq. 1, the PNP equations trivially reduce to the PB equation with the density given by the Boltzmann factor,

(3)

where n₀ denotes a reference density and is the potential energy expressed in a dimensionless form. Using Eq. 3 for n as an integrating factor in Eq. 1, it can be recast into the form

(4)

Under steady-state conditions and assuming a uniform flux J in the z direction, Eq. 4 reduces to 1-D and can be integrated to give

(5)

where the values of  = nexp() at the boundary points z = 0 and L are specified with ₀ and _L, respectively. While Eq. 5 appears to require only the knowledge of the potential in the range [0, L], in fact, there is still a density dependence through Poisson's equation (2). A similar expression for the density can be obtained by integrating Eq. 4 from 0 to z, and using Eq. 5 to eliminate J/D

(6)

Finally, substituting Eq. 6 in Poisson's equation, one obtains an integro-differential equation for the potential in PNP

(7)

This is similar in form to the PB equation, and would reduce to it if _L = ₀, that is, when the electrochemical forces balance out and the system is in equilibrium. In general, there are no known analytical solutions of Eq. 7, and applications of the 1-D PNP to ion channels have to be carried out using numerical methods (Eisenberg, 1996). For future reference, we quote here the trivial solutions of the PNP equations. When the concentration is uniform (n = n₀ everywhere), one simply has Ohm's law

(8)

while in the case of a uniform potential (i.e., no electric forces), the solutions are

(9)

Because a self-consistent analytical solution of PNP equations is not possible, it is natural to look for approximations that will enable such solutions. Even if the potential could be determined in some way, there is still a problem in evaluating the integrals involving its exponential in Eqs. 5-7. In fact, the only known indefinite integral of  exp[f(z)]dz is for f = z, which simply gives back exp(z). This corresponds to the constant field approximation in the Goldman-Hodgkin-Katz theory, and using (z) = ₀ + (_L  0)z/L in Eqs. 5 and 6 yields the following solutions for the flux and density

(10)

(11)

The effect of the electrochemical forces on density, which is not so easy to surmise from Eq. 6, can be seen more clearly from Eq. 11: The density, which varies linearly between the boundary values when there are no electric forces (Eq. 9), exhibits an exponential behavior when there is a uniform field (to be more specific, the density of one type of ions is enhanced while that of the counterions is suppressed relative to the linear case, see Fig. 1 B). Thus the local potential has a significant effect on the cation and anion densities, as in the case of the PB theory. The question is then whether one can calculate this potential correctly in ion channels within the continuum approach using a continuous distribution of charges and mean-field approximation. If we use the BD results in article I as a guide, the answer has to be negative for narrow channels with radius smaller than the Debye length. We have already seen in article I that shielding is largely overestimated in PB theory, and leads to a gross reduction of the potential energy of an ion inside a channel. It is expected that shielding will play a similarly dominant role in the PNP calculations, leading to a largely distorted concentration (and hence current) values in narrow channels when compared to those of the BD simulations.

To test this conjecture, a computer code has been written to solve the PNP equations in three dimensions for a range of channel shapes. In this code, a channel shape is constructed on a rectangular grid, and the PNP equations are solved at the grid points using a finite difference algorithm (see Appendix for details). The input of the program are the channel shape, dielectric constants in the channel and the protein wall, the concentrations and potentials on the reservoir boundaries, the diffusion coefficients of the ions, and the locations and strengths of fixed charges in the channel walls. Once these parameters are specified, the program outputs the concentration and potential throughout the system as well as the ionic currents through the channel. The PNP program is executed on an alpha cluster, where a typical run with 49³ grid points takes 5-10 min. Inclusion of fixed charges in the channel wall roughly doubles the above computation time. When a finer mesh with 99³ grid points is used, the computation time increases by more than an order of magnitude.

Tests of accuracy

View larger version (17K): [in this window] [in a new window]   FIGURE 1   Tests of the accuracy of the PNP code in a cylindrical channel with the same dielectric constant everywhere ( = 80). The length of the channel is 25 Å and results are independent of the radius unless otherwise stated. (A) A comparison of the current flowing through a cylindrical channel of radius 4 Å as found from the PNP code (circles) and from Eq. 8 (solid line). The concentration is set to 300 mM throughout and the potential between the ends of the channel is varied. (B) The concentration profiles of cations and anions in a cylindrical channel when a concentration difference is maintained between the ends (cL = 100 mM and cR = 500 mM). In the absence of an applied potential, the numerical results (diamonds) compare well with the analytical solution from Eq. 9 (dashed line). When a constant field is enforced, the numerical cation (filled circles) and anion (open circles) concentrations again follow closely the analytical results from Eq. 11 (solid lines). The exact PNP concentrations with self-consistent potentials are indicated by the triangles.

Brownian dynamics

Brownian dynamics is relatively new in studies of ion permeation in channels, and not as well known as some other methods. It has been introduced to the field in the review article of Cooper et al. (1985), where a good introduction to the technique in 1-D channel models is given. Despite the encouraging disposition of this article toward BD, it has rarely been taken up in channel studies in the intervening years (Jacobsson and Chiu, 1987; Bek and Jacobsson, 1994). These earlier studies were limited to 1-D channels and could not be expected to model channel-ion interactions correctly. Very recently, we have extended BD simulations to 3-D channels, where all the ion-channel and ion-ion forces are correctly treated according to Poisson's equation (Li et al., 1998; Chung et al., 1998, 1999; Hoyles et al., 1998a). The basic notions of BD simulations are given in the companion article (I). We do not repeat them here but discuss those features of BD specific to ion permeation, which are not mentioned in article I as only the equilibrium situations are considered there.

The Langevin equation (see Eq. 9 in article I) is solved at discrete time steps following the algorithm devised by van Gunsteren and Berendsen (1982). A time step of t = 100 fs is used in the BD simulations in general. In the potassium channel we use a multiple time step algorithm in BD code, as detailed in Chung et al. (1999). A shorter time step of 2 fs is used when an ion is either in the selectivity filter or near the entrance of the channel, where the force acting on an ion changes rapidly. Simulations lasting one or two million time steps are repeated 36 times in the cylindrical channel and 10 to 15 times in the potassium channel to obtain good statistics. A cylindrical reservoir with radius 30 Å is placed at each end of the channel and filled with ions. Its height is adjusted so that the average concentration in the channel/reservoir system remains at the desired value, usually 300 mM represented with a total of 48 ions. This larger value of the concentration than the physiological range (150 mM) is preferred to improve statistics in BD. Initially, the ions in the reservoirs are assigned random positions, with the provision that they do not overlap. Velocities are also assigned randomly according to the Boltzmann distribution. For successive simulations, the final positions and velocities of the ions in the previous simulation are used as initial positions and velocities in the next trial. Current is computed from the number of ions that pass through an imaginary plane at the middle of the channel during a simulation period. To maintain the specified concentrations in the reservoirs, a stochastic boundary is applied: when an ion crosses the channel, say from left to right, an ion of the same species is transplanted from the right reservoir to the left. For this purpose, the ion on the furthermost right-hand side is chosen, and it is placed far left-hand side of the left reservoir, making sure that it does not overlap with another ion. The stochastic boundary trigger points, located at either pore entrance, are checked at each time step of the simulation. Concentration is determined from the time average of the number of ions in a given region. The BD program is executed on a Fujitsu VPP-300 supercomputer. With 48 ions in the system, the typical run time for a simulation period of 1.0 µs (10 million time steps) is ~16 h.

A list of the parameters used in the BD simulations was given in article I. Here we include the mass, diffusion coefficient, and ion radii for K, which was not used in article I:


The same values of the dielectric constants and diffusion coefficients are used in the PNP calculations.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

Comparisons of PNP theory with BD simulations are carried out in cylindrical channels with varying radius, and in a more realistic but complicated model of the potassium channel. The cylindrical channels are the most common geometry used in applications of the PNP approach, and therefore they are used in the majority of tests in the following. Further comparisons are carried out in a model potassium channel that is constructed from its recently revealed structure (Doyle et al., 1998). Unless otherwise stated, the average concentration in the system is kept at 300 mM in both PNP and BD. We note that the Debye length for a 300 mM solution is 5.6 Å.

Cylindrical channels

The cylindrical channel and reservoir system used in both the PNP calculations and the BD simulations is shown in Fig. 2 A. An identical system is used in article I, and as explained there, rounding of the corners is due to the difficulty of solving Poisson's equation with sharp boundaries. The channel radius r is systematically increased from 3 Å to 14 Å, or from 0.5 to 2.5 times the Debye length. The dielectric constants are normally set to 80 in the electrolyte and 2 in the protein wall, except in a control study where _protein = 80 is used to simulate bulk conditions. We use the term "passive channel" to distinguish this non-interacting case from a real channel with _protein = 2. In these comparisons, bare channels (i.e., no fixed charges) are considered first, and then channels with fixed charges in the protein wall.

View larger version (17K): [in this window] [in a new window]   FIGURE 2   (A) Cylindrical channel models used in comparisons of PNP theory with BD simulations. A 3-D channel model is generated by rotating the curve shown about the central axis by 180°. The cylindrical section is 25 Å in length, and the rounded corners have a radius of 5 Å. The radius of the cylinder r is varied from 3 to 14 Å. The reservoir height h is adjusted so as to keep the total (reservoir and channel) volume constant when the radius is changed. (B) The potential energy profiles in a cylindrical channel of radius r = 4 Å when an electric field of 107 V/m is applied in the direction. The dashed and solid lines correspond to the channels with and without fixed charges, respectively. The profile of a passive channel (protein = 80) is indicated by the dotted line.

To ensure that comparisons are carried out in nearly identical situations, we need to match the boundary conditions in PNP with those in BD. Due to the dynamic nature of simulations in BD, there are no unique procedures to implement these conditions. We use a relatively simple strategy here, which will be justified in the control studies below. The applied potential in BD is represented with a uniform electric field (usually E = 10⁷ V/m). The potential difference between the top and bottom boundaries is determined from the potential energy profile of a single ion in the presence of this electric field (see Fig. 2 B), which should yield a reasonable average value when the other ions are present. This potential difference is then used in the PNP calculations. Similarly, the average concentrations in the reservoirs, determined from the number of ions in BD, are implemented in the PNP calculations. There is a slight complication here arising from the fact that concentrations in PNP are specified along the reservoir boundary which, in general, will not match their average values in the reservoir. This happens because any potential drop across the system produces a charge separation. Thus, one needs to find the correct boundary value that will reproduce the desired average concentration in the reservoir. To this end, we first run the PNP program in the absence of any electrolyte to find the potential drop V across the reservoir due to the applied voltage and channel shape. Away from the channel mouth, the charge separation would be expected to balance this potential drop according to the Nernst equation. Thus the ratio of concentrations at the two ends of the reservoir is given by

(12)

which, together with the average concentration, c_av = (c₁ + c₂)/2, determines the appropriate boundary value for concentration. This procedure works well in most cases except when there are fixed charges in the channel or asymmetric solutions are used. These cause further distortions in concentration values that are not taken into account in the above method, with the result that the average concentration in the reservoir does not coincide with the desired value of c_av. In such cases, the PNP runs are iterated until we find the values of c₁ and c₂ that satisfy the Nernst equation (Eq. 12) and the average concentration in the reservoir is equal to the desired value of c_av.

Potential profiles

Control studies

View larger version (29K): [in this window] [in a new window]   FIGURE 3   Comparison of PNP with BD in a passive channel with a radius 14 Å and a symmetric solution of 300 mM in the reservoirs. The ions are driven across the channel with an applied field of E = 107 V/m as indicated in the inset. (A) Concentration profiles of sodium ions (chloride ions exhibit a similar profile, hence not shown here). The channel is divided into 16 layers as shown by the dotted lines in the inset, and each reservoir into 2 layers. The average concentration values in layers are represented by the histograms in BD (reservoir values are shaded) and by the filled circles in PNP. (B) Conductance of Na+ and Cl ions in channels of different radii normalized by the cross-sectional area. The conductance found from the BD simulations is indicated by the filled (sodium) and open (chloride) circles, while those from the PNP theory are shown by the solid lines.

Bare channels

View larger version (36K): [in this window] [in a new window]   FIGURE 4   Comparison of concentration profiles in PNP and BD as in Fig. 3 A but for a real channel (protein = 2) with a radius r = 4 Å. PNP concentrations are shown with filled circles and BD results with the histograms for sodium (A) and chloride (B) ions. A symmetric solution of 300 mM is used and 105 mV is applied between the boundaries.

View larger version (27K): [in this window] [in a new window]   FIGURE 5   Comparison of concentrations in an r = 4 Å channel as in Fig. 4 but with asymmetric solutions (cL = 100 mM and cR = 500 mM) and no applied field.

View larger version (18K): [in this window] [in a new window]   FIGURE 6   Similar to Fig. 4 but shows the changes in the concentration profiles in PNP (solid lines) and BD (dotted lines) as the channel radius is increased progressively from r = 4 to 6, 8, and 12 Å. The concentration of sodium ions is shown in (A) and of chloride ions in (B).

View larger version (15K): [in this window] [in a new window]   FIGURE 7   Normalized conductance of Na+ (A) and Cl (B) ions in a bare channel are plotted against the channel radius as in Fig. 3 B. A symmetric solution of 300 mM and an applied potential of 105 mV are used. The BD results (circles) are fitted by the dotted line and the PNP results (diamonds) by the solid line. Each BD data point is obtained from a 3.6 µs simulation period.

Fixed charges in the channel

View larger version (33K): [in this window] [in a new window]   FIGURE 8   Comparison of concentrations in an r = 4 Å channel as in Fig. 4 but with fixed charges in the protein wall. Eight monopoles, each with charge 0.09e, are distributed around each end of the channel. A symmetric solution of 300 mM and an applied potential of 105 mV are used.

View larger version (23K): [in this window] [in a new window]   FIGURE 9   Concentration profiles in PNP (solid lines) and BD (dotted lines) in cylindrical channels of differing radii as in Fig. 6 but with fixed charges.

View larger version (15K): [in this window] [in a new window]   FIGURE 10   Normalized conductance of Na+ (A) and Cl (B) ions are plotted against the channel radius as in Fig. 7 but for a channel with fixed charges. The BD results (circles), representing a 3.6 µs simulation period, are fitted by the dotted line and the PNP results (diamonds) by the solid line.

View larger version (12K): [in this window] [in a new window]   FIGURE 11   Comparison of I-V curves in PNP (diamonds fitted with solid lines) and BD (circles fitted with dotted lines) in an r = 4 Å channel with fixed charges. An asymmetric solution with cL = 100 mM and cR = 500 mM is used. Each BD point represents a 1 µs simulation period.

View larger version (15K): [in this window] [in a new window]   FIGURE 12   Conductance-concentration curves in PNP (diamonds fitted with solid lines) and BD (circles fitted with dotted lines) in an r = 4 Å channel with fixed charges. Symmetric solutions and an applied potential of 105 mV are used. Each BD data point is derived from a 1 µs simulation period.

View larger version (15K): [in this window] [in a new window]   FIGURE 13   Effect of changing the dielectric constant inside an r = 3 Å channel on sodium and chloride currents. The dielectric constant is kept at 2 in the protein. Both channels without (A) and with (B) fixed charges are considered. A symmetric solution of 300 mM and an applied field of 107 V/m are used. The PNP results are indicated by diamonds fitted with the solid lines. The BD results are shown with circles and mostly vanish.

Potassium channel

The cylindrical channels used in the last section provide only a schematic model for channels. It is of interest to repeat the tests of PNP and BD using a more realistic channel model. For this purpose we use the potassium channel whose crystal structure has been revealed in a recent x-ray study (Doyle et al., 1998). A thorough investigation of this channel using BD is given in Chung et al., 1999. Here we give a minimal description of the model channel necessary for the ensuing discussions. The shape of the channel is shown in Fig. 14 A. A cylindrical reservoir with radius 30 Å and variable length is connected to each end of the channel. The dielectric constants are _water = 80, _protein = 2 as in the cylindrical channels. As shown below, for this narrow channel to conduct it must have fixed charges in the protein wall. These charge groups are modeled as sets of dipoles with fourfold symmetry about the z axis as follows: 1) 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. 2 B) is placed 1 Å from the boundary, the positive pole 1.2 Å away from the negative pole, with their orientation perpendicular to the z axis; 2) 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 Å; 3) four "mouth" dipoles (filled