Test of Poisson-Nernst-Planck theory in ion channels.

The recent Perspectives on Ion Permeation have brought the debate about the applicability of the Poisson-Nernst-Planck (PNP) theory in ion channels to a sharp focus. Despite the differences in opinion, all sides of the debate agree that the mean field approximation in PNP theory needs to be checked by comparison with a more accurate theory; e.g., Brownian dynamics (BD). Clearly, for such a test to be meaningful it has to be carried out in a three-dimensional (3-D) channel. We have been performing 3-D BD simulations in ion channels for the last few years (Li et al. 1998; Chung et al. 1998, Chung et al. 1999; Hoyles et al. 1998), and have been aware of the differences between the two theories. Therefore, we would like to contribute to the debate by providing a simple test of the PNP theory in a cylindrical channel, which appears to be the most common geometry used in applications of PNP.

Due to space limitations, we will not elaborate on either theory, but refer to Chung et al. 1998, Chung et al. 1999 for details of 3-D BD simulations, and Kurnikova et al. 1999 for 3-D–PNP calculations. Reviews of the 1-D BD and PNP can be found, respectively, in Cooper et al. 1985 and Eisenberg 1996. We have written a code similar to the one in Kurnikova et al. 1999 for solving the PNP equations in 3-D. As a control study, the PNP and BD calculations are compared in bulk conditions, and are found to yield the same results for flux and concentration within the computational errors.

A cross section of the channel shape used is shown in Fig. 1 (top). The rounding of corners is required due to the difficulty of solving Poisson's equation with sharp corners. A reservoir with radius 30 Å and variable length is added on both ends of the channel. The length is adjusted so as to keep the concentration fixed at 300 mM when the channel radius is varied. In BD simulations, this concentration is represented by 12 Na+ and 12 Cl− ions in each reservoir. The reason for using a larger value than the physiological range (150 mM) is entirely statistical. Otherwise, almost identical results are obtained for conductance and concentration at 150 mM, once they are normalized to 300 mM. The applied potential in BD is represented with a uniform electric field of E = 107 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. This potential difference is then implemented in the PNP calculations.

Figure 1

Cross section of the cylindrical channel used in the comparisons (top). The length of the channel is 35 Å, including the rounded corners, which have a radius of 5 Å. The dielectric constants are 80 in the channel and 2 in the protein. An electric field of E = 107 V/m is applied as indicated by the arrow. The dotted lines show the layers used in the concentration plot below. (A) Concentration profiles of Na+ ions in a 4-Å radius channel as predicted from PNP theory (•) and BD simulations (histograms). The shaded areas indicate the reservoir values. (B) Conductance of ions, normalized by the cross-sectional area, are plotted against the channel radius. The results obtained from the BD simulations (•, Na+; ○, Cl−) are fitted by the dotted lines. The PNP results are shown by the solid lines.

The potential energy profile of a single ion across the channel exhibits a large barrier due to the repulsive forces from the induced surface charges. This barrier is large enough (≈4 kT) for a 4-Å radius channel to prevent an ion from traversing the channel. When other ions are present in the system, shielding effects might play a role in lowering this barrier and making the channel easier to traverse. The importance of shielding in ion permeation has been amply emphasized in applications of PNP (Eisenberg 1996). However, shielding is a concept developed in continuum theories of bulk electrolytes (e.g., Debye-Hückel theory), and the validity of the mean field approximation inherent in continuum theories is firmly anchored on a comparison of the Debye length of the electrolyte with the system size. A nearly complete screening of an ion's charge occurs at a distance of several Debye lengths in a bulk electrolyte. Noting that the Debye lengths are 7.9 and 5.6 Å for 150 and 300 mM solutions and the channel radii are typically 3–4 Å, one has to worry about the large shielding effects observed in PNP calculations. Are these genuine effects or a chimera caused by an inappropriate application of the PNP theory outside its domain of validity?

To address these questions, we compare in Fig. 1 A the concentration profiles predicted by PNP and BD for Na+ ions in a 4-Å radius channel (a nearly identical result is obtained for Cl− ions). Other than a slight asymmetry in the PNP results caused by the applied potential, the average Na concentration remains around the reservoir value of 300 mM throughout the channel. Since a similar result is obtained for the Cl concentration, we see that almost perfect shielding remains the modus operandi of PNP even though the channel radius is smaller than the Debye length. Perfect shielding implies that there are no induced boundary charges, thus ion-channel interaction is completely ignored in PNP and charge is transported across the channel as if the dielectric boundary did not exist (i.e., ∈protein = 80). Perhaps we should emphasize that there is nothing surprising about the PNP results—these are exactly what one would expect from a theory where ions are represented as a continuous charge density. The question is, if we keep the integrity of ions and take the time average of their motions, would we obtain the same concentrations as predicted by PNP? This question can be answered unambiguously via BD simulations. In Fig. 1 A, we show the BD results (histograms) obtained by averaging the number of Na+ ions in a given layer over the simulation time. Here we have a completely different picture. The Na+ ion concentration drops exponentially as one moves into the channel, and it is more than an order of magnitude smaller than the reservoir values at the middle of the channel. The exponential drop results from the action of the barrier seen by an ion attempting to enter the channel. At times, ions making forays into the channel have sufficient energy (due to fluctuations) to probe the channel interior, but the sharply rising barrier combined with the Boltzmann factor makes this increasingly less probable as the ion gets closer to the center. A similar result is obtained for the Cl− ion concentration. We emphasize that ions enter the channel singly most of the time and not in cation–anion pairs. Thus, the answer to the questions posed above is negative; that is, there is no shielding inside the narrow channels and the representation of ions as continuous charge densities in such situations leads to erroneous results. With increasing channel radius, the discrepancies between the PNP and BD concentrations decrease, and an agreement is achieved at a radius of r ≈ 16 Å. At such a large radius, the effect of the boundary forces on ions becomes insignificant, while the channel is often occupied by counter ions leading to appreciable shielding. Thus the situation becomes similar to the bulk conditions, and once they prevail, the two theories agree as they should.

Since the potential and concentration are determined self-consistently in PNP, the errors committed in concentrations are expected to affect the potential results, and these in turn will lead to inaccuracies in the flux results. To illustrate the magnitude of these errors and how they change with increasing channel size, we plot in Fig. 1 B the normalized conductances obtained from PNP and BD against the radius. The conductance is normalized by dividing it with the cross-sectional area of the cylinder to factor out the trivial increase in flux with the area. The PNP results exhibit a slight reduction with increasing radius, which is due to the access resistance decreasing as 1/r. The BD results, in contrast, start with a zero conductance at r = 3 Å and gradually rise to the PNP levels with increasing radius. The convergence between the two theories occurs at ∼14 Å, which is more than two Debye lengths. The discrepancies observed in conductance results correlate closely with those in concentration discussed above, and thus reaffirm the intimate relationship between the two quantities.

For reasons of clarity, we have discussed above the simplest possible case of a bare channel with symmetric solutions and an applied potential to drive the ions across. In the case of asymmetric solutions, essentially the same conclusions are obtained. Including fixed negative charges in the protein wall, however, introduces an ion-channel interaction that was lacking before and goes some way in reducing the discrepancy between the PNP and BD results. The agreement in concentration is much improved for cations, but not for anions, for which the discrepancy remains about an order of magnitude. The reason is simply that the fixed charges largely cancel the potential barrier on a cation but increase the barrier for an anion. This explains why an artificially small diffusion coefficient for Cl needs to be used in PNP to fit the data. Another problem arises from the ability of the channel to hold arbitrarily large concentrations in PNP: conductance increases linearly with concentration, whereas in BD (as in experiments) it saturates. The BD simulations indicate that the schematic channels considered here are mostly single-ion channels. As pointed out by Levitt 1999, more problems will appear in PNP in the case of multi-ion channels; e.g., potassium channels.