INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

During the last decades continuum theories of electrolytes have found a new niche in the description of physical processes in the salty waters of cells (Weiss, 1996; Evans and Wennerström, 1999). Continuum theories were originally developed for bulk electrolytes early in the century, and their validity has been firmly established since then (Bockris and Reddy, 1970). The more recent applications in biology usually involve mesoscopic systems, and it is not clear from the outset that the assumptions made for bulk solutions are justified for solutions confined to small volumes. Of these, the mean-field approximation that assumes the potential can be determined from a continuous distribution of the mobile charges in an electrolyte is most suspect. The basic question is whether the predicted concentrations in continuum theories, which represent the space average of ion densities, are in accordance with those obtained from Brownian dynamics by time-averaging motions of individual ions. In this respect, the Debye length provides a useful guide. If the system size is much larger than the Debye length, as in the case of large proteins and membrane surfaces, then the mean-field approximation inherent in the continuum theories should be relatively safe. However, membrane pores that transport ions across a cell usually have radii smaller than the Debye length (Hille, 1992), and the use of continuum theories in such systems is questionable. Applications of continuum theories to membrane channels have nevertheless flourished in recent years (for reviews see, for example, Levitt, 1986; Cooper et al., 1985; Eisenberg, 1996, 1999). One reason for the popularity of the continuum theories that deal with concentrations of ions, rather than individual ions, used to be that the alternative methods were computationally intractable. This is still true for molecular dynamics simulations, where motions of both ions and water molecules in the system are traced using Newton's law (Roux and Karplus, 1994). To study ion permeation across a membrane channel using molecular dynamics, one needs supercomputers that are several orders of magnitude faster than currently available. In comparison, the situation with Brownian dynamics (BD), where only the motion of ions are traced, is much better. BD simulations of conductance in realistic three-dimensional geometries are now routinely performed on supercomputers (Li et al., 1998; Chung et al., 1998, 1999; Hoyles et al., 1998a). As stressed in a recent series of commentaries on ion permeation (Levitt, 1999; McClesky, 1999; Miller, 1999; Nonner et al., 1999), the time is ripe for a realistic assessment of continuum theories as models of ion channels, and if they fail the tests, to move on to more accurate theories.

In this and the accompanying article (Corry et al., 2000) we try to provide such a test for two prominent continuum theories: PoissonBoltzmann (PB) in this article and PoissonNernstPlanck in the next one. PB theory has become an important tool for studying proteins and membranes, leading to many insights on the key role played by electrostatic interactions (Honig and Nicholls, 1995). The availability of efficient computer programs for solving the PB equation (Davis and McCammon, 1990; Sharp and Honig, 1990) has increased its use tremendously during the last decade. In ion channels, the PB equation was initially used to include the effects of ionic atmosphere on the potential energy profile of an ion in schematic channel models (Levitt, 1985; Jordan et al., 1989; Cai and Jordan, 1990). More recently, the PB calculation of potential energy profiles has been extended to realistic channel structures in numerous articles (Sankararamakrishnan et al., 1996; Weetman et al., 1997; Adcock et al., 1998; Cheng et al., 1998; Rostovtseva et al., 1998; Sansom et al., 1998; Smejtek et al., 1999; Dieckmann et al., 1999; Ranatunga et al., 1999). In PB calculations, ionic shielding greatly reduces the potential energy of an ion in a channel compared with that of a single ion calculated from Poisson's equation. To assess the reliability of the PB calculations, it is important to check that this shielding effect is not an artefact of the mean-field assumption in a relatively narrow channel environment.

The total electrostatic force acting on an ion inside and near the vicinity of a channel determines its dynamic behavior. Therefore, it is the most important quantity to check in judging the accuracy of the PB theory. Here we test the validity of the mean-field approximation in the PB theory by comparing its predictions for the force on a test ion and potential energy, and concentration profiles with those obtained from BD simulations. BD is eminently suitable for this task because the motion of all the ions in the system are traced individually according to the Langevin equation. Therefore, a long-time average of physical quantities should accurately reflect the actual physical behavior of the system. The main point of this article is demonstrated by using a spherical geometry that serves as a generic example of an electrolyte confined in a small volume. Cylindrical channels with varying radii provide testing grounds for schematic channel models, while a catenary shape similar to that of the acetylcholine receptor channel is used for tests in a more realistic geometry.

	THEORETICAL METHODS

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

PB theory

PB theory provides a classical electrostatic description of a system in which fixed external charges, represented by a density _ex, are surrounded by mobile ions in a dielectric medium. The main assumption of the theory is that at equilibrium, the distribution of mobile ions in the system can be approximated by a continuous charge density, _el, given by the Boltzmann factor

(1)

where n₀ is the bulk (or reference) number density of ions of species  and ze is their charge. Here n₀ (in SI units) is related to concentration c₀ (in moles/liter) by n₀ = 1000N_Ac₀, where N_A is Avogadro's number. The average electric potential (r) in Eq. 1 is obtained from the solution of Poisson's equation

(2)

Combining Eqs. 1 and 2 for a 1:1 electrolyte, which is our main interest here, we obtain the following PB equation:

(3)

Apart from a few special cases this equation cannot be solved analytically. Therefore, a linearized form proposed by Debye and Hückel (1923) has been commonly used in practical applications of the PB equation. Expanding the sinh term in Eq. 3 and keeping only the leading term in  yields the linear PB equation for a bulk electrolyte with no fixed charges (_ex = 0)

(4)

where 1/ is the Debye screening length given by

(5)

At room temperature (T = 298 K) in water ( = 80), the Debye length is related to concentration as ¹ = 3.07/ Å. Although the approximation in Eq. 4 is no longer necessary with the availability of high-speed computers, the intuitive picture of shielding provided by the DebyeHückel theory still plays a useful role. Here we use it to indicate where and why the PB theory may break down. The solution of Eq. 4 in bulk is well known (e.g., McQuarrie, 1976), and yields the following screened Coulomb potential around a central ion of radius a/2

(6)

The radial density of the screening charge p(r) is proportional to this potential

(7)

which is seen to peak at r = 1/ and then decay exponentially. The volume integral of this shielding charge is of interest, and for a sphere of radius r, it is given by

(8)

Equation 8 shows that q(r)/e increases monotonically with r, leading to a 25% screening of the central charge at ~r = 1/, rising to 80% at r = 3/. Thus, for a c₀ = 150 mM electrolyte under bulk conditions, length scales of ~25-30 Å are required for nearly complete screening of an ionic charge. When a boundary is imposed at a smaller distance, the system tries to maintain equilibrium by increasing the counterion concentration in the volume between the ion and the boundary. However, because of the physical size of ions and electrostatic repulsion effects, there is a limit to this increase, and one anticipates that as the ion gets closer to the boundary, the counterion density will eventually diverge from the Boltzmann factor, leading to a much smaller shielding than expected from PB theory. This prediction can be tested directly by comparing the PB results with those obtained from BD simulations, where all ions are treated on an equal footing as particles with a finite size and charge, rather than as a continuous charge density.

For this purpose, we have solved the PB equation (3) numerically by using a finite difference algorithm for various boundaries (see Appendix for details). From the numerical solution of the PB equation, one obtains the potential at discrete grid points. These potential values are then fed into the Boltzmann factor (Eq. 1) to determine the concentration of ions. The components of the force on a test ion at a particular grid point are calculated by using numerical differentiation, from the difference of the potential at two opposing neighboring points in the x, y, and z directions. The PB program is executed on an  cluster, where a typical run with 1 Å grid size takes 5-20 min, depending on the boundary conditions used.

Tests of accuracy

6

Brownian dynamics

Use of Brownian dynamics in studies of ion channels is relatively new. An introduction to the technique for one-dimensional channels is given by Cooper et al. (1985). BD simulations have been extended to realistic three-dimensional channel geometries quite recently (Li et al., 1998; Chung et al., 1998, 1999; Hoyles et al., 1998a). We refer to these articles for further details, and give only a brief description of the method here.

In BD, the motion of individual ions are simulated using the Langevin equation:

(9)

where m_i, q_i, and v_i are the mass, charge and velocity of the ith ion. In Eq. 9, the effect of the surrounding water molecules is represented by an average frictional force with a friction coefficient m_ii, and a stochastic force F_R arising from random collisions. The last term in Eq. 9 is the total electric force acting on the ion due to other ions, fixed and induced surface charges at the channel boundary, and the applied membrane potential. It is computed by solving Poisson's equation for a given channel boundary using an iterative numerical method as detailed in Hoyles et al. (1996, 1998b). Rather than solving Poisson's equation at each time step, which would be computationally prohibitive, a system of lookup tables is used (Hoyles et al., 1998a). The electric field and potential due to one- and two-ion configurations are precalculated at a number of grid points and stored in a set of tables. During simulations, the potential and field at desired points are reconstructed by interpolating between the table entries and using the superposition principle. For this purpose, the total electric potential _i experienced by an ion i is broken into four pieces

(10)

where the sum over j runs over all the other ions in the system. In Eq. 10, _X,i is the external potential due to the applied field, fixed charges in the protein wall, and charges induced by these; _S,i is the self potential due to the surface charges induced by the ion i on the channel boundary; _I,ij is the image potential due to the charges induced by the ion j; and _C,ij is the Coulomb potential due to the ion j. The first three potential terms in Eq. 10 are stored in, respectively, three-, two-, and five-dimensional tables (dimension is reduced by one in the latter two cases by exploiting the azimuthal symmetry of the system's geometry). Similar tables are constructed for each component of the electric field, which are calculated from the gradient of the potential at the grid points.

The Coulomb interaction between two ions is modified by the addition of a repulsive 1/r⁹ potential, which arises from the overlap of their electron clouds (Pauling, 1942)

(11)

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, which is estimated from the ST2 water model used in molecular dynamics as F₀ = 2 ×10¹⁰ N (Stillinger and Rahman, 1974). As demonstrated below, the 1/r⁹ potential emulates the hard-sphere collisions in the primitive model quite well, and hence is adequate for the purpose of comparing PB and BD approaches. However, because the ion pair potential for Na⁺-Cl has a minimum at contact, it leads to some anomalous results in narrow channels (see below). While this is not directly relevant to the main theme of the paper, it is still of interest to point out the source of this anomaly and show that it can be resolved when one uses realistic ionion potentials obtained from molecular dynamics simulations. The hydration forces between two ions add further structure to the ion pair potential in the form of damped oscillations (Guàrdia et al., 1991a, b), which can be approximately represented by

(12)

Here the oscillation length a_w = 2.76 Å is given by the water diameter and the other parameters are determined by fitting Eq. 12 to the potentials of mean force given by Guàrdia et al. (1991a, b). For anion-cation pairs, R_c = r₁ + r₂, but for like ions the contact distance is pushed further to R_c = r₁ + r₂ + 1.6 Å. The origin of the hydration force R is slightly shifted from R_c; by +0.2 Å for like ions and by 0.2 Å otherwise. The exponential drop parameter is determined as a_e = 1 Å for all ion pairs. Finally, the overall strength of the potential is U₀ = 8.5, 2.5 and 1.4 kT for NaCl, NaNa, and ClCl pairs, respectively.

To prevent ions from entering the protein or leaving the system, a hard-wall potential is activated when the ions are at 1 Å distance from the channel or reservoir boundaries, which elastically scatters them. This simple interaction that treats ions as billiard balls with a finite radius is adequate for the purposes of this study. Note that to be precise, the range of the wall potential should correspond to the radius of the ion in question, but for simplicity we have used the same range for all ions. Because only one type of ion (Na⁺) is present in narrow channels in BD simulations, this approximation is not likely to influence the results. A more pertinent question here is the effect of finite ion size in comparisons of BD and PB results. Because ions are represented with a continuous charge density in PB theory, they could occupy the channel right up to the boundary, whereas in BD they are prevented from coming closer than 1 Å to the boundary (that is assuming they carry a point charge at their center). This issue of consistency between the two theories will be addressed in the next section when we compare them in cylindrical channels.

The Langevin equation (9) is solved at discrete time steps following the algorithm devised by van Gunsteren and Berendsen (1982). Throughout, a time step of t = 100 fs is used. Initially, the ions are assigned random positions in the reservoirs, except for the test ion, which is held in a fixed position. 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. A single ion is held in place for a period of 20,000 time steps, while the system reaches equilibrium. After this, the system is allowed to evolve for a further 200,000 or 1,000,000 time steps. At each time step, the force acting on the fixed ion (and other ions) is calculated, and from the time average of these, a value for the force is computed. Each simulation is repeated from between 5 to 16 times to obtain a value for the average force on an ion at each position along the central axis of the system. The duration of simulations is varied from 150 to 300 ns according to the statistical accuracy of the results. The potential profile of an ion is constructed by integrating the force curve along a given path. Average values of the concentration of ions at different points in the system are also obtained for each individual run. The BD program is written in FORTRAN, vectorized and executed on a supercomputer (Fujitsu VPP-300). With 48 ions in the system, the CPU time needed to complete a simulation period of 1.0 µs (10 million time steps) is ~16 h.

A list of the parameters used in the BD simulations is given below:
Dielectric constants:  _water = 80,  _protein = 2; Masses (in kg):  m_Na = 3.8 ×10²⁶,  m_Cl = 5.9 ×10²⁶; Diffusion coefficients (in m² s¹):   D_Na = 1.33 ×10⁹,  D_Cl = 2.03 ×10⁹,   (Note that D is related to the friction constant via D = kT/m); Ion radii (in Å):  r_Na = 0.95,  r_Cl = 1.81; Temperature:  T = 298 K.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

Comparisons of the PB theory with BD simulations are carried out for three different geometries including a sphere, cylindrical channels with varying radius, and a catenary-shaped channel. The cylindrical channels are used in the majority of comparisons because they provide a prototype channel model that have been used in numerous applications of the continuum theories to ion channels. The sphere is included for control studies and pedagogical reasons, and the catenary channel to show the robustness of the results for a more realistic channel shape. Each case is discussed in a separate subsection in the following. We note for future reference that the Debye lengths for 150, 300, and 500 mM solutions are, respectively, 7.9, 5.6, and 4.3 Å.

Electrolyte in a sphere

While our main concern is cylindrical pores, spherical geometry is useful for purposes of control studies in a bulk-like environment, as well as in illustrating the effect of a confining dielectric boundary on shielding of ions in a simple situation. In the following comparisons, a sphere of radius 20 Å containing an electrolyte of concentration c₀ = 500 mM is used. In BD simulations, this concentration is represented by 10 anions and 10 cations, including the test ion. (The systems used here and in the following channel models are always chosen to be electroneutral.) The above choice for concentration is dictated by the BD considerations of having a sufficient number of ions in the system to obtain good statistics, but not too many so as to encumber the simulations. Because the main variable is the distance of the test ion from the boundary, the choice of radius does not have much influence on the results. In order to compare the results with the analytic solutions of the linear PB equation (Eqs. 6 and 7), both the cation and anion radii are taken as 1 Å in the sphere studies. A dielectric constant of  = 80 is used everywhere in bulk simulations. When emulating a protein boundary,  = 2 is used outside the sphere.

In solving the PB equation, we use a sharp spherical boundary around the test ion, which emulates a hard wall potential that prevents its overlap with other ions. Such an infinite potential is not practical to implement in BD simulations; therefore, a 1/r⁹ potential is used instead (Eq. 11), which is both easier to handle and more physical. As seen in Fig. 1 A, the two potentials differ near the contact region and overlap once the ions are slightly separated. As a result of the softer potential used in BD, the ions (especially counterions) are expected to be more broadly distributed near the contact region of the test ion. This is exemplified in Fig. 1 B, where we compare the radial distribution functions g(r) in PB and BD for a bulk electrolyte. Here the PB results are obtained by fixing a test cation at the origin and those of BD by averaging over the ion-ion distributions. To avoid the finite size effects in BD simulations, ion pairs are included in the average only when at least one of them is inside an imaginary r = 10 Å sphere. The linear PB results (not shown) for anion-cation distribution is somewhat higher than the nonlinear one at the maximum but this appears to be mostly due to the finite mesh size used in numerical solutions of PB equation. Otherwise, there is little difference between the linear and nonlinear PB results, especially at larger radii. The broadening of the sharp peak at contact in BD simulations is expected to influence the results at short distances <3 Å. The shifting of counter charge density to smaller radii means that the BD simulations should provide a better shielding at short distances compared to the PB theory. At larger distances, the radial distribution functions overlap, and as far as the force on the test ion is concerned, one should obtain similar results within the two approaches. We note that using a hard-wall potential in BD would have led to larger forces on the test ion at short distances due to less shielding. However, as will be seen in the comparisons below, this issue is mostly irrelevant because the force results in BD closely follow that of a single ion. That is, there is little shielding due to counterions, and therefore details of their interaction with the test ion at short range cannot have much influence on the results.

View larger version (15K): [in this window] [in a new window]   FIGURE 1   (A) Comparison of the hard-wall (solid line) and 1/r9 (dotted line) ion-ion potentials for a positive (U++) and negative (U+) test ion around a fixed positive ion in a 500 mM bulk electrolyte. (B) the resulting radial distribution functions in PB (solid lines) and BD (dotted lines) for anions (g+) and cations (g++) around the fixed ion.

The shortcomings of the continuum theories of electrolytes in confined volumes are most succinctly illustrated in a spherical geometry because it involves a single parameter: the distance of a test ion from the boundary. In Fig. 2 A, the force on a test cation held fixed at a given position is plotted as a function of the radial distance. A single ion (no electrolyte) experiences a repulsive force due to induced surface charges at the sphere boundary. This force is shown by the dashed line for reference purposes. As the ion moves from the center of the sphere toward the dielectric wall, the repulsive force acting on it is seen to increase steeply. The PB calculations (solid line) exhibit the expected results from ionic shielding: the force on the test ion due to the boundary charges is significantly reduced compared to that of a single ion. In contrast, little shielding is observed in BD calculations of the force (filled circles with error bars), which follows quite closely the dashed line for the force on a single ion. The discrepancy between the PB and BD results become appreciable at 8 Å from the boundary, which corresponds to ~2 Debye lengths. As the ion gets closer to the boundary, this discrepancy grows, and at the closest BD simulation point (4 Å), it becomes a factor of 3. We note that even when  = 80 is used outside the sphere, there is a net force on the ion because the presence of a boundary results in an asymmetric distribution of counterions in the radial direction. In both PB and BD, this force is much smaller than the  = 2 case. For example, the force acting on an ion located at 4 Å from the boundary when  = 2 is determined to be 2.5 pN from the PB calculations and 6.8 pN from the BD calculations. The corresponding values when  = 80 are 0.5 and 1.6 pN. Thus the force on the ion is mostly due to the reaction field from the dielectric boundary.

View larger version (15K): [in this window] [in a new window]   FIGURE 2   Test of the PB theory for a 500 mM electrolyte in a sphere of radius 20 Å. (A) Force acting on a fixed cation is plotted against its radial position. The solid curve shows the force obtained from the PB theory and filled circles with error bars show the BD simulation results. The dashed line indicates the force in the case of a single ion (c0 = 0). The error bar on BD data points is one standard error of means and is not shown when it is smaller than the size of the data point. (B) Concentration of mobile ions in a region of constant solid angle (30°) between a fixed cation and the spherical boundary. The average concentration obtained from the PB theory (solid lines) and the BD simulations [open (Cl) and filled (Na+) circles fitted with the dotted lines] are plotted against the radial position of the fixed cation.

The source of the discrepancy is to be sought in the inability of the counterions in the BD simulations to provide the level of shielding observed in the PB theory. To see this more clearly, in Fig. 2 B we compare the anion and cation concentrations in the two theories. The average concentrations in the region between the fixed ion and the boundary are plotted against the radial position of the ion. This region is defined by the conical section between the ion and the boundary with the ion radius as its central axis and subtends a constant solid angle of 30°. Thus as the ion gets closer to the boundary, its volume gets smaller. Concentrations in PB are obtained from the space average of charges in the defined region, whereas in BD they are obtained from the time-average of ions in this region. The PB results are shown with the solid lines, and the BD results are indicated by the open (anions) and filled circles (cations) that are fitted with the dotted lines. As the test ion approaches the boundary, the PB theory predicts a rapid rise in the anion concentration, which is necessitated by the decreasing available volume between the ion and the boundary. The opposite behavior is observed in the BD simulations; that is, the anion concentration actually decreases as the ion gets closer to the boundary. The discrepancy between the predictions of the two theories again becomes appreciable when the ion is ~2 Debye lengths from the boundary. Thus this example explicitly demonstrates when the concentrations in the PB theory start to disagree with the average ion densities obtained from the BD simulations, signaling the break down of the mean-field approximation.

The above results give a clear indication of the operating range of PB theory for an electrolyte confined within a dielectric boundary. While the Boltzmann factor (Eq. 1) puts a limit to the increase in anion concentration (otherwise there would be a perfect shielding with very large concentrations to provide it), it clearly does not capture the whole physical picture. Reflecting on these results, it is clear that the continuum description that distributes the ionic charges over the whole volume is ultimately responsible for the failure of the PB theory. When the integrity of the ionic charges is kept as in the BD simulations, there is an enormous repulsive force on a counterion (due to induced surface charges) as it attempts to enter the narrow region between the test ion and the boundary. This force largely prevents the anions from entering the narrow region, and is responsible for the drop in anion concentration in BD. In PB calculations at 500 mM, an average cell with a grid size of 1 Å contains 1/3000 of a unit charge. Distribution of charge into such small units cuts down the effectiveness of the repulsive force, and hence allows relatively large anion concentrations to occur in the narrow region. While the total negative charge in this region is only ~30% of a unit charge, the PB results indicate that even this small amount could still provide a very effective shielding. This happens because the surface charges induced on the boundary are proportional to 1/r², and therefore those charge elements nearer the boundary can induce proportionately more negative charges on the surface, which cancel the positive charges induced by the test ion more efficiently.

Cylindrical channels

We next consider cylindrical channels with rounded corners. The rounding is necessitated by the fact that sharp corners cause difficulties in numerical solutions of Poisson's equation, and in any case seems to be closer to reality. The dimensions of the channel are outlined in Fig. 3, with the channel obtained by rotating the curve shown in the figure around the symmetry axis. The radius of the channel is varied from 3 Å to 13 Å in the comparisons. The height h of the reservoir is adjusted to keep the volume fixed when the radius is varied. For r = 3 Å, a height of h = 25 Å is used. The dielectric constants are 2 for protein and 80 for water unless otherwise specified. An average concentration of 300 mM is used in all the PB calculations, which is determined from the total cation (or anion) charge in the system as in the case of the sphere. The BD simulations are carried out with a total of 24 Na⁺ and 24 Cl ions, corresponding to an average concentration of 300 mM. The reason for using this higher value instead of the more typical 150 mM is entirely statistical: twice as many ions leads to better accuracy in the BD simulations. The results are hardly sensitive to concentration in BD, and exhibit only a logarithmic dependence in PB calculations. Thus, essentially similar results would be obtained using a concentration of 150 mM.

View larger version (20K): [in this window] [in a new window]   FIGURE 3   Cylindrical channel models used in comparisons of PB theory with BD simulations. A three-dimensional channel model is generated by rotating the cross-section 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 11 Å. The reservoir height h is adjusted so as to keep the total (reservoir and channel) volume constant when the radius is changed.

From the view of the dynamics of an ion in a channel environment, the quantity that is of most interest is the force acting on it at various positions in the channel. In Fig. 4 we compare the PB and BD calculations of the z-component of the force on a test ion as it is moved along the channel axis (only the positive side is shown since the curves are symmetric around z = 0). In BD, a test ion is held at a fixed position on the channel axis, and the z-component of the force acting on it is tabulated at every 10 time steps, which are averaged at the end of the simulation. The ion is then moved to another position along the channel axis, and the measurement is repeated. The shielding effect in PB is seen to lead to a drastic reduction in force compared to the BD result in the r = 3 Å channel (Fig. 4 A). As the channel size is increased, the discrepancy decreases but it remains severalfold (Fig. 4, B and C). Finally, in the r = 11 Å channel, when the force itself becomes quite small, the complete shielding observed in PB theory is reproduced in BD (Fig. 4 D).

View larger version (20K): [in this window] [in a new window]   FIGURE 4   Test of the PB theory for the cylindrical channel shown in Fig. 3. The z-component of the force acting on a fixed cation at various positions along the channel axis is calculated using the PB theory (solid line) and the BD simulations (filled circles fitted with the dotted line). The open circles (where shown) indicate the BD results obtained using the realistic ion-ion potentials given in Eq. 12. The radius of the channel is (A) 3 Å, (B) 4 Å, (C) 7 Å, and (D) 11 Å. The height of the reservoirs is adjusted to keep the concentration fixed at 300 mM in all cases. The force on a single ion (c0 = 0, dashed line) is also shown, for reference purposes.

An interesting observation in Fig. 4, A and B is that the BD results follow rather closely the force on a single ion when the test ion is near the mouth region, indicating that there is absolutely no shielding there, but deviate from it when the ion is further inside the channel. Intuitively, one would have expected the opposite behavior, that is, more shielding of the force as the ion moves out of the channel. To explain the origin of this shielding in BD simulations, we have studied the ion distributions in the r = 3 Å channel when the test ion is at z = 7.5 Å. The channel is found to be devoid of charges except at ~z 11 Å, where there is a net negative unit charge corresponding to an anion trapped in that location. This anion is firmly attached to the test cation, forming a dipole, and thus neutralizing the charge on the test ion to some degree. As pointed out in the Methods section, the short-range ion-ion potential used in the BD simulations has a minimum at contact (see Eq. 11 and Fig. 1 A), and thus is responsible for this anomalous behavior. When realistic ion-pair potentials obtained from molecular dynamics simulations are used instead (Eq. 12), this anomaly disappears and the force on the test ion follows closely that of a single ion, as indicated by open circles in Fig. 4, A and B. Thus, if we discount this anomaly, a fair conclusion of the above study is that shielding does not play any role in the narrow channels with r = 3-4 Å.

Another issue that needs to be addressed in comparisons is the effect of the ion-wall potential (or finite ion size), which is implemented in BD simulations but ignored in PB calculations. This issue of consistency between the two theories and its influence on the results presented can be addressed in two ways. One method is to implement the finite size of ions in PB equation by multiplying the right-hand side of Eq. 3 by a space-dependent function that will exclude the ions from the volume within 1 Å of the boundary (Roux, 1997). The second method is to do the opposite, that is, shrink the activation distance of the hard-wall potential in BD from 1 Å to zero, thus allowing ion centers to come near the boundary. Because in almost all applications of the PB theory to ion channels such finite size effects are not considered and our main purpose is to provide tests for these applications, we prefer to use the second method here. In Fig. 5 we plot the BD results for the force on a cation in an r = 3 Å channel as in Fig. 4 A, but with the activation distance of the hard-wall potential reduced from 1 Å (circles) to 0.5 Å (squares) and 0.1 Å (triangles). It is seen that there are no discernible differences among the various results, with all falling on the force curve obtained from the solution of Poisson's equation for a single ion (dashed line). Thus, even if we ignore the finite size of ions and allow them to access the whole channel volume as in PB theory, they decline to take advantage of the extra space offered. Obviously, the steep increase in image forces as an ion approaches the dielectric boundary makes these regions rather inhospitable places, a fact that is missed by the PB theory because smearing of charges dilutes the effects of the boundary forces. Since the range parameter of the hard-wall potential does not have any influence on the results, we will keep using the more realistic 1 Å range in the rest of the comparisons.

View larger version (14K): [in this window] [in a new window]   FIGURE 5   Effect of changing the activation distance of the hard wall potential in BD simulations. The force on a cation in a r = 3 Å channel is plotted as in Fig. 4 A but for three distance parameters, 1 Å (circles) to 0.5 Å (squares) and 0.1 Å (triangles). The ion-ion interaction from Eq. 12 is used in the simulations. All the BD results follow the single ion results shown by the dashed line. The error bars are not shown to avoid cluttering.

A quantity that can be more directly related to ion permeation is the potential energy profile of an ion which is obtained by integrating the force curves in Fig. 4. We compare the PB and BD profiles in Fig. 6 for r = 3, 4, and 7 Å channels (the r = 11 Å results are not shown because both profiles are too small on the scale of the graphs). In the PB case, shielding reduces the energy barrier seen by a single ion by roughly an order of magnitude, virtually obliterating it. Ionic atmosphere is seen to provide some shielding in BD by lowering the energy barrier, though in the case of the narrow channels (r = 3-4 Å), this is caused by the short-range potential used as explained above. If one uses a realistic ion-ion potential, the energy barrier goes back to that of a single ion (open circles in Fig. 6, A and B). In the r = 7 Å channel, shielding is seen to reduce the barrier of a single ion by more than half. Since a similar result is obtained with the realistic ion-ion potential, this is likely to be a genuine result. Nevertheless, the barrier in BD remains larger than the PB result, pointing to a sizable discrepancy despite the reduction in the potential energy values.

View larger version (17K): [in this window] [in a new window]   FIGURE 6   The potential energy profiles in PB (solid line) and BD (dotted line) obtained by integrating the force curves in Fig. 4. The dashed line shows the profile of a single ion. The open circles (where shown) indicate the BD results with the realistic ion-ion potentials (Eq. 12).

Compared to the sphere results, the discrepancy between the two theories is more accentuated in the cylindrical channels because the access of counterions to a narrow cylinder is further hindered in BD, while no such hindrance occurs in PB theory. To quantify this statement, we compare the anion (Fig. 7 A) and cation (Fig. 7 B) concentrations predicted by PB (solid line) and BD (bars) theories for an r = 3 Å channel when the test ion is located near the pore mouth (z = 12.5 Å). In PB calculations, both anions and cations uniformly occupy the channel at about the average concentration, except near the test ion when the former shoots to very large values and the latter dips to zero as expected. This difference in the anion and cation concentrations leads to a net screening charge of 0.61e in the channel. In stark contrast, both anions and cations are completely excluded from the channel interior in BD. While there is some excess of counterions near the channel entrance, these only amount to 0.01e, which is too small to provide any shielding as seen from the force at z = 12.5 Å in Fig. 4 A. We note that a constant anion concentration of 300 mM throughout the channel would correspond to a total charge of 0.21e. Thus the amount of anion charge is increased severalfold compared to the background in PB, while it remains negligibly small in BD.

View larger version (18K): [in this window] [in a new window]   FIGURE 7   Variation of the average concentration along an r = 3 Å channel for cations (A), and anions (B) when a cation is fixed on the z axis at z = 12.5 Å (where the channel starts curving). In BD, the channel is divided into 32 layers and the average value of the concentration is calculated at each layer. The PB concentrations are indicated by the solid curve and the BD ones by the bar graph.

Rather than repeating the above study for each channel size, which is not very informative, we demonstrate the changes in concentration by plotting the total screening charge in the channel as a function of its radius (Fig. 8 A). This study is carried out for a cation fixed at z = 12.5 Å, where the force on an ion is at a maximum. The total screening charge in PB remains nearly constant with the increasing radius, the slight increase being due to approaching bulk conditions (note that the screening charge in the channel remains less than e because the channel volume is limited to z = ±15 Å). In BD, this charge is negligible at r = 3 Å but it steadily rises with r, converging to the PB value at ~r = 11 Å or 2 Debye lengths. As shown in Fig. 8 B, the force on the test ion at z = 12.5 Å correlates very well with the screening charge results in (A). The force in BD initially coincides with that of a single ion at r = 3 Å (no shielding), and with increasing radius, it gradually converges to the PB values at around 2 Debye lengths. This study establishes the domain of validity of PB theory for channels as 2 Debye lengths, below which the underlying mean-field approximation breaks down to an increasingly larger degree with decreasing radius.

View larger version (14K): [in this window] [in a new window]   FIGURE 8   Pore size dependence of the screening charge and force on a cation held at z = 12.5 Å. (A) The net screening charge in the channel (from z = 15 to 15 Å) is plotted as a function of the channel radius. The PB results are shown by the solid line and the BD values by the filled circles fitted with the dotted line. (B) Force on the cation as the channel radius is increased. The PB (solid line), BD (filled circles fitted with a dotted line) and single ion results (dashed line) are indicated in the figure.

So far we have considered only the central axis in comparisons, which may give the impression that an agreement between the PB and BD results can be obtained in the larger-size channels (Fig. 4 D). However, the central axis is a rather special place where the forces from the boundary charges are at a minimum, and the shielding effects in BD are maximized due to the azimuthal symmetry. From the sphere results it is expected that as the ion is moved toward the boundary, discrepancies between the two theories will resurface. This is quite obvious for the radial component of the force but not so for the z-component. In Fig. 9 we present comparisons of the z-component of the force on a test ion similar to Fig. 4 D (reproduced at the top) but along an axis that is offset from the z axis by 4 Å (Fig. 9 B) and 8 Å (Fig. 9 C). Shielding effects are again overestimated in PB theory compared to BD as the ion approaches the boundary. To see this more clearly, we show in Fig. 10 how the z-component of the force changes as the ion is moved radially from the center to the boundary. Up to 3 Å from the center the two theories agree, that is, shielding of the force in PB theory is reproduced in BD. But after that, shielding progressively weakens in BD in contrast to PB theory, which provides a very good shielding right up to the boundary. Thus, even in large channels the predictions of the PB theory are bound to fail as one approaches the channel walls. Such discrepancies in large channels are relevant, for example, in calculation of concentrations near a binding site, but not in ion transport, as ions tend to stay near the channel axis where the radial force is minimum (Li et al., 1998; Corry et al., 2000).

View larger version (19K): [in this window] [in a new window]   FIGURE 9   Comparison of the z-component of the force on a test ion in an r = 11 Å channel when it is offset from the central axis by r = 4 Å (middle) and r = 8 Å (bottom). The top figure is the same as Fig. 4 D.

View larger version (16K): [in this window] [in a new window]   FIGURE 10   Comparison of the z-component of the force on a test ion in an r = 11 Å channel as it is moved radially from the center to the channel boundary at z = 8.75 Å.

The BD results so far clearly indicate that narrow channels with radii 3-4 Å are pretty inhospitable places for ions regardless of their background concentration. Therefore, for ion permeation to take place it is essential to reduce the energy barrier of a bare channel by placing fixed charges of opposite sign on the protein wall. To test the PB theory in this more realistic case, we place a set of negative charges in the walls near the each end of an r = 3 Å channel. Eight monopoles with charges 0.09e are spread evenly around the channel circumference at z = 12.5 Å and z = 12.5 Å, where the channel starts curving. The PB, BD, and single ion results for the z-component of the force on a test cation (as in Fig. 4 A) are compared in Fig. 11 A. The fixed negative charges on the channel wall reduce the presence of counterions in the channel and the associated shielding, and hence lead to much larger forces in PB theory compared to Fig. 4 A, in better agreement with the BD results. The fourfold discrepancy observed in the bare channel (Fig. 4 A) is now reduced to about a factor of 2. The fixed charges also prevent the occurrence of anomalous shielding in the channel interior; the BD results now closely track those of the single ion, reinforcing the earlier conclusion that no shielding is possible inside a narrow channel. As the test ion is moved away from the channel, shielding becomes more and more effective and the force in BD goes gradually from the single ion curve toward the PB result. In Fig. 11 B we show the potential energy profiles obtained from the force curves in Fig. 11 A. Besides the usual discrepancy between PB and BD theories, perhaps a paradoxical result is that shielding actually increases the energy barrier in BD compared to that of a single ion. The reason for this ironic result can be seen from Fig. 11 A; shielding operates when the ion is outside the channel where the force is attractive, but not inside when it is repulsive.

View larger version (17K): [in this window] [in a new window]   FIGURE 11   Effect of placing fixed charges in the channel wall in an r = 3 Å channel. (A) Force on a cation as in Fig. 4 A but with fixed charges. (B) Potential energy profiles as in Fig. 6 A but with fixed charges.

As a final study in cylindrical channels, we consider the possibility that the dielectric constant inside the channel may be smaller than 80, especially in narrow channels. This can be implemented in a straightforward manner in the PB algorithm where a three-dimensional grid is used, but it is not so easy in the BD simulations where a boundary element method is used in solving Poisson's equation. This problem has been tackled in previous BD simulations (Chung et al., 1998, 1999) by using the reduced value of the dielectric constant, _c, in both the channel and reservoir, and including the neglected Born energy difference between the channel-reservoir configurations with _c-80 and _c-_c as a short-range energy barrier at the channel entrances. We refer to the above references for details of this implementation in the BD program. The Born energy difference for a 3 Å channel is calculated using the PB program at zero concentration, which gives a barrier height of 3.5 kT for _c = 40 and 11.8 kT for _c = 20. In PB calculations the change in  is implemented in five equal steps from the channel entrance at z = 17.5 Å to z = 12.5 Å, and similarly at the other end. The potential energy profiles for _c = 20 and 40 are compared in Fig. 12. The barrier height for a single ion increases roughly as 1/_c, and a similar trend is seen in BD. At lower _c, BD results deviate more from those of single ion because of the appearance of shielding at the mouth region. We attribute this to the stronger Coulomb attraction between the test ion and counterions, which increases as 1/_c. The corresponding increase in barrier height is much faster in PB theory, so that the discrepancy with BD gets smaller with decreasing _c (but stays severalfold in any case). This faster increase of barrier height in PB is related to the loss of shielding inside the channel with reduction in , which affects the PB results but not BD. Nevertheless, the overall conclusion remains the same as before; greater shielding in PB results in much lower energy barriers compared to BD.

View larger version (17K): [in this window] [in a new window]   FIGURE 12   Effect of changing the dielectric constant in the channel, c, on the potential profile of a test ion in an r = 3 Å channel without fixed charges. (A) c = 20; (B) c = 40.

Catenary channel

The above study in cylindrical channels gives a good idea about the expected working range of the PB theory. To demonstrate the robustness of those conclusions, we repeat the force and potential energy calculations in a more realistic channel geometry with vestibules. This catenary-shaped channel is generated by rotating the closed curve shown in Fig. 13 about the axis of symmetry. The vestibule of this channel is similar in shape to that visible in the electron microscope pictures of the acetylcholine receptor channel (Toyoshima and Unwin, 1988), making this a better approximation of a real biological ion channel. The vestibules are generated by a hyperbolic cosine function, z = a cosh(x/a), where a = 4.87 Å. The entrance to the vestibule has a fixed radius of 13 Å. Two such identical vestibules are connected to a cylindrical transmembrane segment of radius 4 Å and length 10 Å. It is assumed for convenience that the vestibules have the same shape and size, although the electron microscope images show the extracellular vestibule to be larger than the intracellular vestibule.

View larger version (19K): [in this window] [in a new window]   FIGURE 13   Diagram showing the cross-section of the catenary geometry that approximates the shape of the acetylcholine receptor channel. A three-dimensional channel is generated by rotating the curves about the central axis by 180°. Vestibules at each side of the membrane are constructed using a hyperbolic cosine function, y = a cosh(x/a) where a = 4.87 Å. The radius at the entrance of the vestibule is 13 Å and at the cylindrical transmembrane segment 4 Å. Cylindrical reservoirs (not shown), 30 Å in radius and 22 Å in height, are attached to the vestibules.

In Fig. 14 A we show the z-component of the force as the test ion is moved along the central axis of the channel. The concentration is maintained at 300 mM in both the PB calculations and the BD simulations. As before, PB calculations are shown by the solid line, the BD results are indicated by the filled circles which are fitted by the dotted line, and the dashed line shows the force on a single ion. The BD calculations of force closely track the single ion results in the narrow parts of the channel (up to z = 10 Å), and reinforce the earlier conclusion on impossibility of shielding in narrow parts of the channel. There is a large discrepancy between the PB and BD results in this region as in Fig. 4, A and B. As the ion is moved further along the z axis, the channel expands and shielding becomes more and more effective in BD. This is reflected in the force values in BD gradually moving from the single ion curve to the PB results in the z = 10-30 Å range. The potential energy profiles obtained from the force curves in Fig. 14 A are shown in Fig. 14 B. Shielding is seen to have reduced the energy barrier of a single ion by 40% in BD; however, the barrier in BD is still three times larger than the PB result. Thus, in a channel with vestibules, shielding definitely plays some role but its effect is nowhere near the PB predictions, where shielding demolishes the barrier presented to a single ion for all practical purposes.

View larger version (18K): [in this window] [in a new window]   FIGURE 14   (A) The z-component of the force on a cation for a 300 mM electrolyte in a catenary channel is plotted against its axial position. The force obtained from the PB theory is shown by the solid curve and the BD results are indicated by filled circles with error bars fitted with the dotted line. The dashed line indicates the force in the case of a single ion (c0 = 0). (B) The potential energy profiles obtained from the force curves in (A).

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION THEORETICAL METHODS RESULTS AND DISCUSSION CONCLUSIONS APPENDIX REFERENCES

The comparisons of PB theory with BD simulations in various configurations provide clear answers on the range of validity of the former. When the distance of an ion from the channel wall is less than 1 Debye length, the PB calculations largely overestimate the shielding effects and cannot be expected to give reliable values of the force on and potential energy of an ion. The convergence of the PB and BD results occurs when the ion's distance from the channel wall is ~2 Debye lengths, depending on the quantity and the geometry considered. Because the radii of membrane channels are typically smaller than the Debye length, the PB theory cannot be used to obtain reliable estimates of electrostatic force and potential energy of an ion in the channel environment.

Our BD results demonstrate that if the radial profile of a channel is less than the Debye length throughout, it is unlikely to contain any counterions. This conclusion is especially reinforced in realistic channel configurations where fixed charges of opposite sign, which are necessary for ion permeation, make it virtually impossible for any counterion to enter the channel (cf. Fig. 11). For such channels, it is clearly better to use Poisson's equation rather than PB, as no shielding due to ionic atmosphere is possible. This conclusion appears ironical in the historical context of the field because PB theory was advanced as an improvement of Poisson's equation in ion channels. Channels whose radial profiles exhibit large variations are more difficult to reconcile with the existing continuum theories because each is valid in a limited range. For such channels, BD simulations certainly offer a more reliable method for calculations of forces and potentials. Nevertheless, if one insists on using a continuum description, one could presumably extrapolate from Poisson's to the PB equation as the channel widens by using the BD results as a guide.