This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Corry, B. Articles by Chung, S.-H. Search for Related Content PubMed PubMed Citation Articles by Corry, B. Articles by Chung, S.-H.

Dielectric Self-Energy in Poisson-Boltzmann and Poisson-Nernst-Planck Models of Ion Channels

Ben Corry ^*, Serdar Kuyucak ^* and Shin-Ho Chung 

^* Department of Theoretical Physics, Research School of Physical Sciences and Department of Physics, Australian National University, Canberra, Australia

Correspondence: Address reprint requests to Ben Corry, Dept. of Theoretical Physics, Research School of Physical Sciences, Australian National University, Canberra, A.C.T. 0200, Australia. Tel.: 6-126-125-2052; Fax: 6-126-247-2792; E-mail: ben.corry{at}anu.edu.au.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
We demonstrated previously that the two continuum theories widely used in modeling biological ion channels give unreliable results when the radius of the conduit is less than two Debye lengths. The reason for this failure is the neglect of surface charges on the protein wall induced by permeating ions. Here we attempt to improve the accuracy of the Poisson-Boltzmann and Poisson-Nernst-Planck theories, when applied to channel-like environments, by including a specific dielectric self-energy term to overcome spurious shielding effects inherent in these theories. By comparing results with Brownian dynamics simulations, we show that the inclusion of an additional term in the equations yields significant qualitative improvements. The modified theories perform well in very wide and very narrow channels, but are less successful at intermediate sizes. The situation is worse in multi-ion channels because of the inability of the continuum theories to handle the ion-to-ion interactions correctly. Thus, further work is required if these continuum theories are to be reliably salvaged for quantitative studies of biological ion channels in all situations.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Continuum theories of electrolytes, such as the Poisson-Boltzmann (PB) and Poisson-Nernst-Planck (PNP) equations, have been widely used in studies of ion channels during the last two decades. For example, the PB equation has been employed in calculation of the forces and energy profiles seen by ions in channels (e.g., Levitt, 1985; Jordan et al., 1989; Weetman et al., 1997; Adcock et al., 1998; Cheng et al., 1998), and to determine the ionization states of side chains (Roux et al., 2000; Ranatunga et al., 2001) under equilibrium conditions. Similarly, the PNP equations have been used to investigate the permeation properties of many channels (for review see Eisenberg, 1999; and for recent applications see Kurnikova et al., 1999; Cardenas et al., 2000; Nonner et al., 2000; Hollerbach et al., 2000; Hollerbach and Eisenberg, 2002). Recently the application of continuum theories to narrow ion channels (i.e., radius smaller than the Debye length of the electrolyte) was called into question, as they were found to overestimate electrolytic shielding effects (Moy et al., 2000, Corry et al., 2000, Graf et al., 2000). The PB equation predicts that high concentrations of counterions accumulate around a fixed test ion in a narrow channel, which act to eliminate the repulsive image forces due to the polarization charges induced at the channel walls. Such high concentrations of counterions around a test ion are not seen when the mobile ions are treated explicitly as in Brownian dynamics (BD) simulations. In BD, the image forces felt by individual counterions prevent them from easily entering the channel and thereby suppress the large degree of shielding seen in the PB theory. Similarly, in PNP calculations the currents passing through narrow channels are found to be much larger than in corresponding BD simulations, again as a direct consequence of higher ionic concentrations found in the former method.

Extensive comparisons with BD simulations indicate that the problems in both continuum theories arise from the fact that the dielectric self-energy contribution to the ions' potential energy is not properly taken into account (Moy et al., 2000, Corry et al., 2000). Because the self-energy is proportional to the square of the charge, it is positive for all ions, and it would repel them from the dielectric channel boundaries leading to a suppression of their concentrations. In numerical solutions of the PB and PNP equations, however, the charge of mobile ions is distributed throughout the channel and only a fraction of an ionic charge is found in a given region. This leads to a reduced amount of induced charges on the channel boundary, and the quadratic dependence of the self-energy on charge means that its effect is significantly underestimated. In reality, an ion carrying a unit charge of e induces charges on the channel boundary in an all-or-nothing way—either a large image force is induced (repelling the ion away from the boundary) or there is no ion and thus no force. The lack of dielectric self-energy was accentuated in many of the PNP calculations carried out previously, in particular those with no fixed charges in the channel walls (Corry et al., 2000). In these PNP results equal amounts of positive and negative electrolyte were found in the channel, and so there was no net charge at any given point. Thus, no charge was induced along the dielectric boundaries. This was in contrast to the BD simulations in which the positive and negative ions entered the channel at different times, and so each induced a large amount of surface charges and felt a large repulsive force from the channel boundaries. As shown by Corry et al., 2000, the failure to treat induced charges correctly in standard PNP models leads to currents much larger than those observed in BD.

Mean field theories of electrolytes have been commonly used to determine charge transport in semiconductors (Haensch, 1991; Assad et al., 1998) as well as to describe the electric field around globular proteins (Davis and McCammon, 1990; Sharp and Honig, 1990). Thus it may be surprising that the lack of dielectric self-energy has only been noted as a problem relatively recently. However, generally the dimensions of these systems are much greater than the Debye length. In such cases the effects of image forces can be shielded out by the surrounding electrolyte. Image forces are also more significant inside channels where the ion is surrounded by dielectric boundaries rather than facing them on just one side (Kuyucak et al., 2001). Thus, there has been little need to deal with the interaction with dielectric boundaries more carefully in semiconductors or around globular proteins as is required in ion channels.

Given the computational simplicity of the continuum theories compared to simulation approaches and their widespread use, it would be desirable to find a solution to the problem of dielectric self-energy highlighted above. A possible approach is to include an explicit self-energy term in the formalism that would mimic the effect of the induced surface charges on ions entering the channel. Such a term would prevent the buildup of large counterion concentrations around a test ion and may eradicate the spurious shielding effects seen in the continuum models. This approach to correct the continuum theories is being pursued elsewhere (Schuss et al., 2001; Mamonov et al., 2002; Coalson et al., 2002), but no explicit tests of its performance have been published. Here we examine the consequences of including a dielectric self-energy term in the PB and PNP formalisms. The forces calculated using the modified Poisson-Boltzmann and the currents found from the modified Poisson-Nernst-Planck equations are again compared with Brownian dynamics simulations to assess whether the suggested modification improves the performance of the continuum approaches in ion channels.

Another potential problem in applications of the continuum theories to ion channels is the spurious self-energy that arises from continuous distribution of an ion's charge. Distributing an ionic charge across a region leads to interactions of the charge with itself that are not present when the charge is localized at a single point. In single-ion channels, this is a relatively minor problem compared to that of dielectric self-energy, and perhaps because of this, it has not been considered in earlier literature. The significance of this effect, however, grows rapidly with the number of ions resident in a channel. Therefore, to understand and interpret the results obtained in multi-ion channels we will also discuss this spurious self-energy contribution in continuum theories in some detail.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Modified PB equation
In Poisson-Boltzmann theory, the mobile ions are treated as a continuous charge density distributed according to the Boltzmann factor

(1)

Here n₀ is a background (or reference) density of ions of species , U is the electrostatic potential energy of an ion with charge ze, and (r) is the average potential at the position r determined from the solution of Poisson's equation

(2)

where _ex represents all the other external charge sources in the system such as fixed charges in the protein. The PB equation follows upon substituting the density from Eq. 1 into Eq. 2.

A crucial assumption of PB theory is that at a given point, the same average potential acts on all ions regardless of their valence, and hence the potential energy of anions and cations have opposite signs. This is a reasonable assumption to make in a homogeneous system but it cannot be justified in an inhomogeneous one with dielectric boundaries. Near a water-protein interface, a mobile ion with charge ze induces polarization charges of the same sign on the boundary. These charges, in turn, generate a reaction potential _R that acts to repel the ion from the boundary. The potential energy due to this reaction field, called the dielectric self-energy, is given by

(3)

Since _R is itself proportional to ze, U_R depends on the square of the ionic charge, and hence it is always positive regardless of the valence of ions. Thus the dielectric self-energy contribution to the total potential energy of an ion is the same for cations and anions, which is incompatible with the assumption in Eq. 1. The upshot is that the PB equation cannot lead to a suppression of ionic densities near a dielectric boundary because it fails to take into account the dielectric self-energy contribution. A similar problem occurs in the Gouy-Chapman solutions of the PB equation where counterions build up at a charged electrode. In that case, this unphysical behavior is avoided by introducing a Stern or Helmholtz layer that excludes ions from the immediate neighborhood of the boundary. While the presence of this layer is motivated by the finite size of ions or the effect of the solvent layer, in practice its thickness is used as a free parameter to fit the experimental data. In any case, the physics of repulsive self-energy is very different and cannot be modeled by such a simple exclusion zone for ions.

The above discussion about the lack of dielectric self-energy in the PB equation also suggests an apparently simple remedy: modify the Boltzmann factor in Eq. 1 by adding a dielectric self-energy term to the potential energy

(4)

so that the modified PB equation in a 1:1 electrolyte becomes

(5)

As an ion approaches a dielectric boundary, U_R grows rapidly, and the exponential factor in Eq. 5 provides a natural mechanism to suppress the ion densities. The difficulty in this scheme, of course, lies in the calculation of U_R in a many-body system in a self-consistent manner. For a single ion, the self-energy is well-defined—one just needs to solve Poisson's equation for the single ion and substitute the computed reaction potential in Eq. 3. In a dilute electrolyte solution, where the effect of the reaction field created on one ion by another can be neglected, this simple recipe should be quite adequate for the self-energy correction. However, at higher concentrations, the approximation of the reaction potential with that of a single ion becomes problematic. Ion-to-ion correlations can be included into continuum electrolyte theories, such as has been done to describe the formation and behavior of highly charged ‘macro-ions’ and DNA molecules (e.g., Grosberg et al., 2002; Gelbart et al., 2000). In most cases these studies do not include the effect one ion has on another due to the reaction potential created at a dielectric interface as the calculations are carried out in a region of uniform dielectric constant, although some include dielectric boundaries with simple geometries (Nguyen et al., 2000). In continuum theories, incorporating the effect that the reaction field that is created by one ion has on another with more complicated geometries is a difficult problem. Thus, it is of interest to see whether the inclusion of the self-energy correction at the single-ion level (i.e., only the influence the reaction field created by an ion has on itself) improves the accuracy of the PB solutions sufficiently so that they can be used in ion channels with some confidence. For this purpose, we compare the results of the modified PB equation with those obtained from BD simulations.

The modified PB equation is solved using a finite difference method as described previously for the standard PB equation (Moy et al., 2000). The values of dielectric self-energy are calculated by solving Poisson's equation using the same technique, but with all the electrolyte concentrations set to zero. A unit charge is placed at each of the grid points used in the finite difference computation and the calculated values of U_R are stored on the computer in a table. Entries from this table are utilized during the solution of the modified PB equation (Eq. 5).

Modified PNP equations
The conduction of ions in an electrolyte solution is described by the Nernst-Planck equation

(6)

where D is the diffusion coefficient of the ions of species . When Eq. 6 is combined with Poisson's equation (Eq. 2), one obtains the PNP equations, which provide a consistent framework for the description of ion transport in homogenous systems. As in the case of the PB theory, the presence of dielectric boundaries creates problems for the application of PNP to narrow channels because the dielectric self-energy of ions is not properly taken into account.

The PNP equations can be modified in a similar manner to the PB equation by including a specific dielectric self-energy term in the potential energy of an ion:

(7)

The analogy between Eq. 7 and the modified PB equation can be seen most clearly when the current vanishes. Then Eq. 7 can be easily integrated yielding the modified Boltzmann factor given in Eq. 4. Here we solve Eq. 7 simultaneously with Eq. 2 using a finite difference procedure as described previously for the PNP equations (Corry et al., 2000). The values of dielectric self-energy are implemented in this process using a table as in the solution of the modified PB equation.

Self-energy in continuum theories
In continuum electrostatics, self-energy refers to the intrinsic potential energy of a charge (or charge distribution) with itself. This is different from the dielectric self-energy discussed above, which arises from the interaction of a charge with a dielectric boundary. Because an ion's self-energy is constant and has no effect on the results, it is routinely ignored in electrostatic calculations involving discrete charges, and only the Coulomb interactions among the charges are considered. In continuum theories, however, charges are distributed continuously and it is not possible to isolate the self-energy from the Coulomb interaction. Thus when Poisson's equation is solved in the PB and PNP theories, the calculated potential energy of ions necessarily contains a spurious self-energy contribution. This is not a concern in typical applications of continuum theories, which involve a large number of ions (N >> 1) in a nearly neutral electrolyte solution. First, the self-energy is proportional to N whereas Coulomb interaction scales as N², so its relative contribution would be negligible for large N. Secondly, near electroneutrality in a large system means that the self-energy due to any excess charge has to be small. To give an example, when a charge q is uniformly spread in a water-filled sphere of radius R, the self-energy associated with this distribution is given by

(8)

For a unit charge with a dielectric constant of _w = 80, this yields U_s = 4.2 / R kT where R is in Å. Thus, in a large system (e.g., R > 100 Å), self-energy due to an excess ionic charge is completely negligible.

Ion channels do not satisfy either condition. They contain a rather small volume of electrolyte which is far from being electroneutral—at least for the majority of channels that are either cation- or anion-selective. Because the PB and PNP equations strive to maintain electroneutrality, this was not an apparent problem in earlier applications of the continuum theories to ion channels. The introduction of a dielectric self-energy term, however, leads to a suppression of the counterions in the channel and thereby exposes the problem of self-energy more conspicuously. To appreciate the magnitude of this effect, we show in Fig. 1 its influence in the simple case of electrolyte contained within a sphere. In Fig. 1 A, a single ionic charge is confined inside a 4 Å radius sphere with = 80 everywhere so that there are no dielectric boundaries. A single ion should, thus, move freely within the sphere occupying all space evenly. In contrast, the average concentration found using the (standard) PB equation (Fig. 1 A) exhibits a central depression, which is simply caused by the Coulomb repulsion among the charges distributed at all grid points. The potential energy of this charge distribution is 1 kT. This self-energy is, of course, entirely spurious because the potential energy of an isolated ion is zero.

View larger version (15K): [in this window] [in a new window]   FIGURE 1  Spurious self-energy contribution to the total potential energy when the ionic charges are distributed continuously. In both cases the charge is contained inside a 4-Å radius sphere and = 80 is used everywhere. (A) The concentration found using the standard PB equation is plotted against the radial position when there is a net charge of e in the sphere. (B) The energy of the system as the charge in the sphere is increased. Results are shown as found from the PB equation at 298 K (solid line), 0 K (dash-dot line), and as calculated with discrete ions at 0 K (dashed line).

While channels are cylindrical and do not absolutely confine ions, we expect the order-of-magnitude estimates presented above will have a bearing on the results of the modified PB and PNP equations, particularly in ion channels occupied by more than one ion. Whenever there is a net charge buildup inside the channel, the repulsion arising from the spurious self-energy will tend to push the electrolyte out of the channel, leading to a lowering of the concentration inside. This additional repulsion mechanism is required to understand some of the results presented below, especially in the case of multi-ion channels.

Brownian dynamics
Brownian dynamics simulations are carried out to compare results for the force on ions in the channel and the current across the channel with those found in the modified PB and PNP equations. BD provides a useful comparison as the mobile ions in the electrolyte are treated explicitly, whereas the dielectric boundaries are handled using continuum electrostatics in a similar manner to the continuum approach. Thus differences between BD and continuum theory results arise directly from using discrete ions rather than a continuous charge distribution.

The use of Brownian dynamics in ion channel simulations has been described previously in our previous articles (Chung et al., 1999; Moy et al., 2000; Corry et al., 2000, 2002; Kuyucak et al., 2001), and indeed some of the same data are utilized here as in those earlier studies. Very briefly, in BD the motions of individual ions are simulated using the Langevin equation in which the effects of surrounding water molecules are represented by frictional and random forces. The electric forces acting on the ion at each time step, including Coulomb interactions, interactions with the dielectric channel boundary, and any applied membrane potential are calculated by solving Poisson's equation using an iterative boundary element method (Hoyles et al., 1998). Additional forces that represent the short range interaction of ions, hydration effects, and interactions with the channel boundary are also employed (Corry et al., 2001). Ions are initially given random positions in reservoirs attached to each end of the channel, and the Langevin equation is solved at discrete time steps of 100 fs following a third-order algorithm (van Gunsteren and Berendsen, 1982). In the multi-ion channels, a shorter time step of 2 fs is utilized for ions inside the channel.

When comparing forces on an ion in the channel with the solutions of the modified PB equations, an ion is held at a fixed position in the channel while the other ions are free to move throughout the system. The force acting on the ion is averaged over the simulation period. When comparing with the modified PNP equations, all ions are free to move.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Because our aim is to see how much the inclusion of the dielectric self-energy term improves the continuum theories, we will compare the results of the modified PB and PNP equations with those obtained from BD simulations as well as from the standard PB and PNP equations. As the latter two theories have been compared extensively previously (Moy et al., 2000; Corry et al., 2000), we will rely on some of those results here. Therefore, the same schematic cylindrical channel with rounded corners is adapted in most of the comparisons (Fig. 2). The channel has a fixed length of 35 Å and its radius is varied from r = 3 to 13 Å. Two reservoirs with a radius of 30 Å and a variable height are attached to the channel. The height of the reservoir varies with the radius of the channel. For example, the height is h = 25 Å when the pore radius is 3 Å. For larger channel radii, h is adjusted to keep the entire volume of the system constant. The dielectric constants are set to 2 in the protein and 80 for the water. An average concentration of 300 mM is used in all cases, corresponding to 24 Na⁺ and 24 Cl^- ions in the BD simulations. The Debye length for this concentration is 5.6 Å.

View larger version (21K): [in this window] [in a new window]   FIGURE 2  Cylindrical channel models used in comparisons of PB and PNP 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 curvature of 5 Å, bringing the total length of the channel to 35 Å. The radius of the cylinder r is varied from 3 to 13 Å. The reservoir height h is adjusted so as to keep the total (reservoir and channel) volume constant when the radius is changed.

View larger version (17K): [in this window] [in a new window]   FIGURE 3  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 modified PB (MPB) results are shown by the solid line, standard PB by the dashed line, and the BD values by the solid circles fitted with the dash-dot line. (B) The axial component of the force on the test ion normalized by the force on a single ion (calculated by solving Poisson's equation) is plotted as the channel radius is increased. Symbols are as in A.

Biological ion channels contain fixed charges in the protein walls that create valence selectivity. For completeness, we examine the more realistic case of a r = 4 Å channel with negative charges at the entrances. These negative charges cancel the energy barriers encountered by cations, while roughly doubling their height for anions, thus making the channel cation-selective. Here, we place a ring of eight negative monopoles each with charge -0.09 e spread evenly near the channel mouths at z = ±12.5 Å and set 1 Å inside the boundary. In Fig. 4 A we plot the potential profile through the channel as found from the PB equation (upper dashed line), modified PB equation (MPB, solid line), BD simulations (dash-dot line) and Poisson's equation with no electrolyte (lower dashed line). Here there is no test ion in the channel as this has already been considered above. The potential found from the modified PB equation shows a significant improvement compared to the standard PB results. Interestingly, although the modified PB potential displays a considerable drop, it does not overshoot the BD result as in Fig. 3.

View larger version (22K): [in this window] [in a new window]   FIGURE 4  (A) Potential profiles found in a 4-Å radius cylindrical channel with charges in the protein and no fixed test ion found using the standard PB equation (upper dashed line), modified PB equation (solid line), and BD (dash-dot line) with 300 mM NaCl solution in the baths. The potential found from Poisson's equation with no electrolyte is shown by the lower dashed line. (B) Concentration profiles corresponding to the results in A.

We next consider the gramicidin A channel, which offers a more realistic case study. A model of the channel is constructed from the recent NMR data (Ketchem et al., 1997), including the partial charges of all the protein atoms as described previously (Edwards et al., 2002). In Fig. 5, A and B, we show the potential profile and concentration of cations found in this channel in a similar manner to Fig. 4. In the absence of an electrolyte, the potential has a deep well created by the partial charges in the protein, which act to attract cations into the channel. Using the standard PB equation, the potential is flattened out because, as seen in Fig. 5 B, a large concentration of cations enters the channel and almost cancels the Coulomb potential of the partial charges. When the modified PB equation is used, however, we find that this potential profile remains largely unchanged from the solution of Poisson's equation because the dielectric self-energy term leads to a large suppression of the cation concentration in the channel. The potential found from the modified PB equation is remarkably similar to that found from the BD simulations, the two agreeing within a few percent. Comparing the concentration profiles in Fig. 5 B, we find that in both cases the concentration of ions in the channel is quite low. Thus in the gramicidin A channel, the inclusion of the dielectric self-energy term in the PB equation yields a dramatic improvement. This may not be surprising as the channel is very narrow and is found to be occupied only 10% of the time by a single cation. This means that the dielectric repulsion from the channel boundary is a dominant effect, and it can be approximated quite well by the repulsion felt by a single ion.

View larger version (15K): [in this window] [in a new window]   FIGURE 5  (A) Potential and (B) concentration profiles as in Fig. 4 except in the gramicidin A channel. A concentration of 500 mM KCl is used in both cases. The shape of the channel is shown in the inset.

View larger version (16K): [in this window] [in a new window]   FIGURE 6  (A) Potential and (B) concentration profiles as in Fig. 4 except in the KcsA potassium channel. A concentration of 300 mM KCl is used in both cases. The shape of the channel is indicated in the inset.

Modified PNP equations
As a first test case, we again consider the schematic cylindrical channel shown in Fig. 2 without any charges in the protein. Following the tests carried out in a previous article (Corry et al., 2000), we calculate the current passing through this channel using the modified PNP equations and compare them to currents derived from BD simulations. All calculations are carried out using a symmetrical 300 mM NaCl solution and an applied potential of 105 mV between the reservoir ends. In Fig. 7, we plot the channel conductance for cylindrical channels as a function of the channel radius. Here the conductance is normalized by the cross-sectional area of the channel to factor out the trivial increase in current with increasing area. The old PNP results of Corry et al. (2000) are indicated by the dashed lines and show the slight downward trend created by access resistance effects. The currents in narrow channels found from BD (data points fitted by dash-dot lines) are well below the PNP values due to the repulsive self-energy barrier retarding the permeation of ions through the channel. A similar suppression occurs in modified PNP, when the dielectric self-energy term is included explicitly in the PNP equations. Indeed, the currents carried by both cations and anions are essentially zero in the 3-Å radius channel in agreement with the BD results. The normalized conductance climbs as the channel radius is increased in a similar manner to that seen in BD. There are some differences in intermediate radii channels, but the differences are remarkably small (35%). It appears that the inclusion of an explicit dielectric self-energy term in the PNP equations can reproduce the BD results reasonably well in these bare channels.

View larger version (17K): [in this window] [in a new window]   FIGURE 7  Conductance of Na+ (A) and Cl- (B) ions for bare (no fixed charge) channels of varying radii normalized by the cross-sectional area. The results of the modified PNP equations (solid lines), standard PNP (dashed lines), and BD simulations (data points fitted by dash-dot lines) are shown. The ions are driven across the channel with an applied field of 105 mV between the reservoir ends and a 300 mM NaCl solution is maintained in the reservoirs.

View larger version (17K): [in this window] [in a new window]   FIGURE 8  Normalized conductance in cylindrical channels with fixed charges in the channel walls; otherwise, as in Fig. 7.

View larger version (22K): [in this window] [in a new window]   FIGURE 9  Concentration profiles for (A) Na+ and (B) Cl- ion in a 4-Å radius cylindrical channel with fixed charges as found from the modified PNP equations (solid line), standard PNP (dashed line), and BD simulations (bars). The pore is divided into 16 equal segments along the channel axis, and the average concentration in each segment is plotted.

View larger version (15K): [in this window] [in a new window]   FIGURE 10  Current-concentration relationships for (A) CaCl2 and (B) NaCl in the L-type calcium channel. Results found from BD simulations (data points fitted by dash-dot line), the standard PNP equations (dashed line), and the modified PNP equations (solid line) are shown. The shape of the channel and the locations of two of the four glutamate residues (squares), and the negative ends of the intracellular helix dipoles (diamonds) are shown in the inset. Diffusion coefficients used are 1.33 and 0.79 x 10-9 m2 s-1 for sodium and calcium.

Another difference between the modified PNP and BD results worth noting is the different number of ions that are found in the channel. The average number of cations and anions in the channel with either 150 mM CaCl₂ or 150 mM NaCl that were found using the different theories are shown in Table 1. In BD simulations the conduction of calcium ions is essentially a two-ion process. One calcium ion resides in the selectivity filter, and a second must enter for conduction to take place. Indeed, the most common situation is for two ions to reside in the selectivity filter, and the rate-limiting step in conduction is for one of these ions to exit the selectivity filter. Thus, it is found that there are, on average, almost two ions in the channel. The modified PNP equations, however, predict just over one ion in the channel. It also predicts a more realistic number of sodium ions in the channel, but overall the inclusion of the dielectric self-energy term does not improve the results of the theory when considering the number of ions in the channel.

View larger version (31K): [in this window] [in a new window]   FIGURE 11  Average number of ions in layers comprising the right-half of the calcium channel as found from BD simulations (bars) as well as the standard (dashed line) and modified (solid line) PNP equations using 150 mM CaCl2. The section of the channel for which the concentration is plotted is indicated by the dashed box in the inset.

Finally we examine currents in a potassium channel. Rather than using the full-atom model described earlier when using the Poisson-Boltzmann equation, we have employed a simplified channel that gives similar currents as described elsewhere (Chung et al., 2002). The reason for using this model is that it is simpler to include in the PNP routines and we can also adjust the dimensions of the model more easily, such as, for example, widening the intracellular pore diameter to examine how such alterations vary the results of the PNP theory. The shape of the channel is shown at the top of Fig. 12.

View larger version (18K): [in this window] [in a new window]   FIGURE 12  Current passing through the simplified potassium channel model plotted against the intracellular pore radius of the channel as found from BD simulations (data points fitted by dash-dot line), the standard (dashed line), and modified (solid line) PNP equations. The shape of the channel is shown in the inset.

When we consider the average number of ions in the channel, we find that the modified PNP equations again predict too few ions as shown in Table 2, particularly in the narrower channels. Whereas the BD simulations find 3.5 ions in the narrow 2-Å radius channel, the modified PNP equations find one less: only 2.5 ions. Although this may not seem important if the currents are reasonably accurate, knowing the number of ions involved in the permeation process is essential if we are to understand the process properly. Thus the inclusion of the dielectric self-energy term certainly improves the PNP results, but substantial errors are still encountered in the predicted currents and ion concentrations.

	DISCUSSION AND CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

When using the modified PB equation, the greatest improvement is obtained in the narrowest single-ion channels studied. This is because the dielectric self-energy term has the largest effect in such channels, and its inclusion in the PB equation almost eradicates ion concentrations in the channel in accordance with the BD simulations. Agreement between the modified PB theory and BD simulations is also retained in wide channels with r > 2 Debye lengths. However, at intermediate channel radii the agreement is not as good, and errors of up to 50% are common. The validity of the modified PNP theory in single-ion channels mimics that of the modified PB equation. Accurate results are obtained in very narrow and very wide channels, but errors are present at intermediate radii. Inclusion of the dielectric self-energy term yields a general qualitative improvement, but where the standard continuum equations overestimate electrolyte shielding, the modified continuum equations underestimate it in many cases.

The overall performance of the modified continuum equations is worse in the multi-ion potassium and calcium channels considered here. In the KcsA potassium channel, the concentration obtained from the modified PB equation is found to be too low in the selectivity filter. The qualitative shape and magnitude of the potential also exhibit significant discrepancies when compared to the BD results. The modified PNP equations do better when applied to a highly charged calcium channel model. The predicted calcium currents are a little too high and the sodium currents too low, but, overall the inclusion of dielectric self-energy yields more realistic current values. In both the KcsA and calcium channel models, the channel is highly charged and there are multiple ions resident in the channel. The interactions between the discrete ions are expected to play an important role in ion dynamics, which is difficult to describe in the continuum picture. For example, the magnitude of the dielectric repulsion acting on an ion will be modified in the presence of a second ion. Because the dielectric self-energy correction is introduced at the single-ion level, this effect is likely to be underestimated in a multi-ion channel. In addition, the spurious self-energy associated with a continuous charge distribution becomes more important in a multi-ion channel, although this may somewhat counteract the lack of dielectric self-energy. The incorrect numbers of ions inhabiting the multi-ion channels predicted by the continuum theories may be a manifestation of the inaccuracies involved in describing interactions between ions in the channel. The problems highlighted in these and earlier studies also indicate that the validity of using the PB equation for calculating pK_A values in ion channels needs to be assessed.

More sophisticated statistical mechanical descriptions of electrolytes are being developed, but as far as we are aware, none of these theories attempt to take into account the dielectric self-energy term in ion channels. The primary applications of density functional theory in inhomogeneous systems are concerned with electrolytes near a highly charged plane (Henderson et al., 2000) or hard-sphere fluids in a cavity (Gonzalez et al., 1997), neither of which is directly relevant to ion channels. In the former, the image forces are simply ignored (considering the high charge on the planes, this is presumably a good approximation), whereas in the latter there are no image forces.

Recently, there have been a number of attempts to build finite ion size effects into continuum models of ion channels using the mean spherical approximation or density functional theory (Nonner et al., 2000; Goulding et al., 2000; Boda et al., 2002). It is worth considering whether the discrepancies between the modified continuum theories and BD simulations could be improved by doing this. The BD simulations differ from the continuum models by allowing finite sizes for ions, as well as by including them as discrete charges which induce surface charges. Thus, it is possible that this is creating discrepancies between the results of the models, as well as the lack of dielectric surface charge. However, the inclusion of ion sizes would only worsen the level of agreement seen between the modified PB theory and BD simulations seen here. This is clear if we look at the cases where the theories disagree. As stated above, the modified PB theory tends to overestimate the size of the repulsive forces acting on ions in the channel. This means that the concentrations of ions are below that found from BD. Taking into account the finite size of ions in the continuum theory would act to limit how densely charge can be packed into the channel. This, of course, would tend to decrease the concentration in the channel further and worsen the agreement with BD simulations. Thus, it is important that the effects of dielectric self-energy are considered in future mean field studies of ion channels.

The continuum theories are designed for use in systems containing a large number of ions, which is obviously not satisfied in ion channels. Therefore, their application to ion channels needs to be continuously monitored to check against spurious and nonphysical effects. Here we have shown that attempts to improve the PNP and PB equations by including a specific dielectric self-energy term work reasonably well in some cases, but not all. This appears to be a promising direction for improvement, but further work is required to remove all the shortcomings highlighted in this work. Until a more sophisticated continuum description can be developed that overcomes the problems raised here, simulation methods remain as more reliable tools for investigation of ion channels. Brownian and molecular dynamics are now well-developed for this purpose. They can be used to predict basic permeation properties of channels such as conductance that used to be the forte of the continuum theories. In addition, they allow exploration of phenomena such as ionic mixtures and ion selectivity that cannot be studied using continuum approaches.

	ACKNOWLEDGEMENTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
The calculations upon which this work is based were carried out using the Fujitsu VPP-300, the Linux alpha cluster of the ANU Supercomputer Facility, and the Compaq AlphaServer SC of the Australian Partnership for Advanced Computing.

This work is supported by grants from the Australian Research Council, and the National Health and Medical Research Council of Australia.

Submitted on October 24, 2002; accepted for publication January 31, 2003.

	REFERENCES

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION AND CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Adcock, C., G. R. Smith, and M. S. P. Sansom. 1998. Electrostatics and the ion selectivity of ligand-gated channels. Biophys. J. 75:1211–1222.[Abstract/Free Full Text]

Assad, F., K. Banoo, and M. Lundstrom. 1998. The drift-diffusion equation revisited. Sol. Sta. Electron. 42:283–295.

Boda, B., D. Henderson, and D. D. Busath. 2002. Monte Carlo study of the selectivity of calcium channels: improved geometrical model. Mol. Phys. 100:2361–2368.

Cardenas, A. E., R. D. Coalson, and M. G. Kurnikova. 2000. Three-dimensional Poisson-Nernst-Planck theory studies: influence of membrane electrostatics on gramicidin A channel conductance. Biophys. J. 79:80–93.[Abstract/Free Full Text]

Cheng, W., C. X. Wang, W. Z. Chen, Y. W. Xu, and Y. Y. Shi. 1998. Investigating the dielectric effects of channel pore water on the electrostatic barriers of the permeation ion by the finite difference Poisson-Boltzmann method. Eur. Biophys. J. 27:105–112.[Medline]

Chung, S. H., T. W. Allen, M. Hoyles, and S. Kuyucak. 1999. Permeation of ions across the potassium channel: Brownian dynamics studies. Biophys. J. 77:2517–2533.[Abstract/Free Full Text]

Chung, S. H., T. W. Allen, and S. Kuyucak. 2002. Conducting-state properties of the KcsA potassium channel from molecular and Brownian dynamics simulations. Biophys. J. 82:628–645.[Abstract/Free Full Text]

Coalson, R. D., M. G. Kurnikova, P. Graf, and A. Nitzan. 2002. Brownian dynamics studies of ion permeation through biological channels. Biophys. J. 82:340a.

Corry, B., S. Kuyucak, and S. H. Chung. 2000. Tests of continuum theories as models of ion channels. II. Poisson-Nernst-Planck theory versus Brownian dynamics. Biophys. J. 78:2364–2381.[Abstract/Free Full Text]

Corry, B., T. W. Allen, S. Kuyucak, and S. H. Chung. 2001. Mechanisms of permeation and selectivity in calcium channels. Biophys. J. 80:195–214.[Abstract/Free Full Text]

Corry, B., M. Hoyles, T. W. Allen, M. Walker, S. Kuyucak, and S. H. Chung. 2002. Reservoir boundaries in Brownian dynamics simulations of ion channels. Biophys. J. 82:1975–1984.[Abstract/Free Full Text]

Davis, M. E., and J. A. McCammon. 1990. Electrostatics in biomolecular structure and dynamics. Chem. Rev. 90:509–521.

Doyle, D. A., J. M. Cabral, R. A. Pfuetzner, A. Kuo, J. M. Gulbis, S. L. Cohen, B. T. Chait, and R. MacKinnon. 1998. The structure of the potassium channel: molecular basis of K⁺ conduction and selectivity. Science. 280:69–77.[Abstract/Free Full Text]

Edwards, S., B. Corry, S. Kuyucak, and S. H. Chung. 2002. Continuum electrostatics fails to describe ion permeation in the gramicidin channel. Biophys. J. 83:1348–1360.[Abstract/Free Full Text]

Eisenberg, R. S. 1999. From structure to function in open ionic channels. J. Membr. Biol. 171:1–24.[Medline]

Gelbart, W. M., R. F. Bruinsma, P. A. Pincus, and V. A. Parsegian. 2000. DNA-inspired electrostatics. Phys. Today. 53:38–44.

Gonzalez, A., J. A. White, F. L. Roman, S. Velasco, and R. Evans. 1997. Density-functional theory for small systems: hard spheres in a closed spherical cavity. Phys. Rev. Lett. 79:2466–2469.

Goulding, D., J. P. Hansen, and S. Melchionna. 2000. Size selectivity in narrow pores. Phys. Rev. Lett. 85:1132–1135.[Medline]

Graf, P., A. Nitzan, M. G. Kurnikova, and R. D. Coalson. 2000. A dynamic lattice Monte Carlo model of ion transport in inhomogeneous dielectric environments: method and implementation. J. Phys. Chem. B. 104:12324–12338.

Grosberg, A. Y., T. T. Nguyen, and B. I. Shlovskii. 2002. The physics of charge inversion in chemical and biological systems. Rev. Mod. Phys. 74:329–345.

Haensch, W. 1991. The Drift-Diffusion Equation and its Applications in MOSFET Modeling. Springer-Verlag, New York.

Henderson, D., P. Bryk, S. Sokolowski, and D. T. Wasan. 2000. Density-functional theory of an electrolyte confined by thin charged walls. Phys. Rev. E. 61:3896–3903.

Hollerbach, U., D. P. Chen, D. D. Busath, and B. Eisenberg. 2000. Predicting function from structure using the Poisson-Nernst-Planck equations: Sodium current in the gramicidin A channel. Langmuir. 16:5509–5514.

Hollerbach, U., and R. S. Eisenberg. 2002. Concentration-dependent shielding of electrostatic potentials inside the gramicidin A channels. Langmuir. 18:3626–3631.

Hoyles, M., S. Kuyucak, and S. H. Chung. 1998. Solutions of Poisson's equation in channel-like geometries. Comp. Phys. Commun. 115:45–68.

Jordan, P. C., R. J. Bacquet, J. A. McCammon, and P. Tran. 1989. How electrolyte shielding influences the electrical potential in transmembrane ion channels. Biophys. J. 55:1041–1052.[Abstract/Free Full Text]

Ketchem, R. R., B. Roux, and T. A. Cross. 1997. High-resolution polypeptide structure in a lamellar phase lipid environment from solid state NMR derived orientational constraints. Structure. 5:1655–1669.[Medline]

Kurnikova, M. G., R. D. Coalson, P. Graf, and A. Nitzan. 1999. A lattice relaxation algorithm for three-dimensional Poisson-Nernst-Planck theory with application to ion transport through the gramicidin A channel. Biophys. J. 76:642–656.[Abstract/Free Full Text]

Kuyucak, S., O. S. Andersen, and S. H. Chung. 2001. Models of permeation in ion channels. Rep. Prog. Phys. 64:1427–1472.

Levitt, D. G. 1985. Strong electrolyte continuum theory solution for equilibrium profiles, diffusion limitation, and conductance in charged ion channels. Biophys. J. 48:19–31.[Abstract/Free Full Text]

Mamonov, A., R. D. Coalson, P. Graf, A. Nitzan, and M. G. Kurnikova. 2002. Closer look at the dielectric barrier in the narrow channel, real problem or theoretical artifact? Biophys. J. 82:209a.

Moy, G., B. Corry, S. Kuyucak, and S. H. Chung. 2000. Tests of continuum theories as models of ion channels. I. Poisson-Boltzmann theory versus Brownian dynamics. Biophys. J. 78:2349–2363.[Abstract/Free Full Text]

Nguyen, T. T., A. Y. Grosberg, and B. I. Scklovskii. 2000. Screening of a charged particle by multivalent counterions in salty water: strong charge inversion. J. Chem. Phys. 113:1110–1125.

Nonner, W., L. Catacuzzeno, and B. Eisenberg. 2000. Binding and selectivity in L-type calcium channels: a mean spherical approximation. Biophys. J. 79:1976–1992.[Abstract/Free Full Text]

Ranatunga, K. M., I. H. Shrivastava, G. R. Smith, and M. S. P. Sansom. 2001. Side-chain ionization states in a potassium channel. Biophys. J. 80:1210–1219.[Abstract/Free Full Text]

Roux, B., S. Bernèche, and W. Im. 2000. Ion channels, permeation, and electrostatics: insight into the function of KcsA. Biochemistry. 39:13295–13306.[Medline]

Schuss, Z., B. Nadler, and R. S. Eisenberg. 2001. Derivation of Poisson and Nernst-Planck equations in a bath and channel from a molecular model. Phys. Rev. E. 64:0361161–14.

Sharp, K. A., and B. Honig. 1990. Electrostatic interactions in macromolecules: theory and applications. Annu. Rev. Biophys. Biophys. Chem. 19:301–332.[Medline]

van Gunsteren, W. F., and H. J. C. Berendsen. 1982. Algorithms for Brownian dynamics. Mol. Phys. 45:637–647.

Weetman, P., S. Goldman, and C. G. Gray. 1997. Use of Poisson-Boltzmann equation to estimate the electrostatic free energy barrier for dielectric models of biological ion channels. J. Phys. Chem. 101:6073–6078.

This article has been cited by other articles:

				H. Gong, G. Hocky, and K. F. Freed Influence of nonlinear electrostatics on transfer energies between liquid phases: Charge burial is far less expensive than Born model PNAS, August 12, 2008; 105(32): 11146 - 11151. [Abstract] [Full Text] [PDF]