Coupling Poisson-Nernst-Planck and Density

Functional Theory to Calculate Ion Flux

Wolfgang Nonner  Dirk Gillespie  Robert S. Eisenberg

May 28, 2002

1 Introduction

The movement of ions through open ion channels is one of the many interesting physical chemistry problems presented by living cells. In this paper we describe a system to model such transport, using it as an example of a theory that may be applied to other problems of ion flux.

Cells and cell organelles are enveloped by lipid membranes that are nearly impermeable to physiological ions (mostly Na^{+}, K^{+}, Ca^{2+}, and Cl^{−}). One mechanism for ions to move across these membranes is through open ion channels, proteins embedded in the membrane that form ion-selective pores. This type of charge movement conducts electrical signals down nerves and initiates muscle contraction, to list just two of the many physiological functions of ion channels [1].

Ion channels are proteins with functional groups that often are the side chains of charged amino acids. Enough data on these functional groups and even crystallographic channel structures are now available to apply specific theories of permeation and selectivity. With such information and theories it is possible to reverse engineer the devices that biologists observe, not only with the goal of understanding their natural function, but also of controlling their function beyond biology. In this paper we describe an engineering approach to selective ion conduction through ion channels. Specifically, we develop a one-dimensional model of the movement of ions between two baths of fixed ionic concentrations and applied electrostatic potentials. The baths are connected by a single open ion channel whose functional groups are represented as ions confined to the channel. In this model (unlike some we have studied [2], [3], [4]) ions are charged, hard spheres immersed in a hard-sphere solvent and uniform, continuum dielectric. Particle transport is described as friction-limited drift-diffusion with all excess chemical potentials described by the Density Functional Theory (DFT) of Rosenfeld and colleagues [5], [6], [7], [8], [9]. As an example, we present flux calculations for a simple model of a biological calcium channel. Although our specific goal is to model ion channels, the system is generally applicable to one-dimensional modeling of ion transport.

2 Geometry

I do not understand the slab geometry. I need to see the figure and BE SURE that EVERY reference to slab geometry makes sense. This is VERY important because otherwise it will surely attract the attention if not anger of the reviewers.

The geometry of two baths connected by a single ion channel embedded in a membrane is shown in Fig. 1. This is the geometry of the most direct experiments [10], [11]. One inherent difficulty in modeling ion channels is the different length scales, with the nanometer-scale channel connected to the millimeter-scale baths where boundary conditions essential for the transport are maintained. Because the flux is controlled in the channel, the resolution must be high there. However, the boundary conditions of the problem must be applied far away from the channel to model the problem correctly. This problem is three dimensional, but a reasonable one-dimensional approximation can be developed that still includes the effects of the finite volume of ions and functional groups of the channel protein. In this section we derive an approximation that is fine grained in the flux-limiting region (the channel) and coarse grained in the baths. This approximation allows the problem to be solved on a desktop computer.

Starting in the baths, consider spherical shells that terminate at the lipid and are perpendicular to the lipid as well as the x-axis (Fig. 1). The flux of particles is perpendicular to these surfaces and therefore flux densities scale inversely with the area A (x) of these shells. (The shells are indexed by their intersection with the x-axis.) This spherical geometry is exact for each bath alone: ions traveling to and from the channel must converge to a hemisphere given by the capture radius of the channel. As we move from the baths into the channel, we continue to assume that particles travel perpendicular to these spherical shells. Below we will argue that this is an appropriate approximation inside the channel where the spherical shells are flat, cross-sectional disks. In the transition regions this approximation is also valid. For example, in the calculation of the electric field in these transition regions, use of the spherical shells is equivalent to approximating a flat-bottomed cone by a spherical cone [12].

In the channel itself, it is important to consider the nature of the boundary between the pore (permeation pathway) and the channel protein that surrounds it. If there were a smooth, hard wall at the pore/protein interface that maintained a fixed volume for the pore, then large radial packing effects would result in the pore [13]. However, real proteins are not smooth or rigid, and therefore we take the opposite view of a hard-wall interface: the atoms of the channel protein are another liquid, forming a statistical interface with the contents of the pore and providing a confining normal pressure on the pore contents. Such a boundary is not described by a geometric surface at a specific location, but rather by an average cross-sectional area or pore volume that may be constant or vary with conditions depending on the compressibility of the protein [14]. Here we consider the simple case that the protein maintains a constant average cross-section. Comment: the word ‘variable’ might imply ‘vary with time’ and that would be confusing. I find the word ‘variable’ often causes this problem for readers.

Using this description of the interface, we can approximately describe the output variables (the concentrations and the electrostatic potential) as well as the input variables (especially the dielectric coefficient) in the slab geometry inside the channel (that is, they are constant throughout each cross-sectional disk and on to infinity in the radial direction) because:

1. From the hard-sphere perspective, the protein makes the pore effectively wider than its physical radius. Any radial space demands made by the particles in the pore will affect (radially) not only the contents of the pore, but also the particles of the protein; the protein atoms compete for space with the particles in the pore, as opposed to a hard wall that would merely confine the permeating particles to the pore. Seen as only a hard-sphere liquid, the combination of the pore and protein is approximately one continuous liquid and because it is radially relatively large, it may be approximated in the slab geometry.

2. By describing the electrostatics in the slab geometry, each cross-sectional disk is an equipotential surface and therefore there is no electric flux out of the pore; this is equivalent to the situation that the exterior of the pore has a dielectric coefficient of zero. The electric polarizability of the pore contents and the protein are not known, in particular on the time scale of ion conduction. However, because the lipid has a dielectric coefficient of approximately 2, as we move from the pore through the protein to the lipid, the zero dielectric coefficient approximation becomes better and better. In general, by using the slab geometry we imply that as far as we can describe the polarization around permeating ions by an effective dielectric coefficient, we can approximate the polarization originating from regions away from the pore as negligible. Furthermore we imply that we can neglect the electric field reaching out into the lipid as negligible or even zero.

3. In the cases to which we will apply this model (such as the XXX calcium channel we will consider later), there are such large concentrations of ions inside the channel that the screening length becomes less than half the average radius of the pore. Brownian dynamics studies have shown that under these conditions the slab geometry can be used in mean-field theories like the one we consider here [15].

This representation of the pore/protein interface is not a means to obtain a one-dimensional approximation. In future two- or three-dimensional calculations the protein must be modeled to include the two important features discussed above: the hard-sphere interactions extend through the protein while the electrostatic interactions are essentially confined to the pore. IMPORTANT: I do not understand. If the word ‘above’ means previously, I could not find the discussions. If ‘above’ means later, I also could not find the features discussed. This will surely attract criticism from reviewers. COMMENT: I think the ideas that hard sphere interactions extend through the protein while electrostatic interactions are essentially confined to the protein are VERY important. So there should be paragraphs that say:

"Hard sphere interactions extend through the pore because....."

"Electrostatic interactions are confined to the pore (in contrast to hard sphere interactions) because ..."

3 Modeling ion flux

We describe steady-state flux by the constitutive law

−J_{i} = \frac{1}{kT} D_{i} (x) A (x) ρ_{i} (x) \frac{dμ_{i}}{dx}  (1)

and the continuity equation

\frac{dJ_{i}}{dx} = 0;  (2)

that is, for particle species i, the (constant) flux J_{i} is proportional to the particle density ρ_{i} (x) and the gradient of the electrochemical potential μ_{i} (x). A (x) is the area of the spherical shells described in the previous section. This model was recently derived from a Poisson-Langevin system [16].

The electrochemical potential consists of the ideal component μ_{i}^{id} (x), the excess component μ_{i}^{ex} (x) from particle interactions, and the concentration-independent component μ_{i}^{0} (x):

μ_{i} (x) = μ_{i}^{0} (x) + μ_{i}^{id} (x) + μ_{i}^{ex} (x)  (3)

with

μ_{i}^{id} (x) = z_{i}eϕ (x) + kT ln [\frac{ρ_{i} (x)}{ρ_{scale}}]  (4)

where ϕ (x) is the local electrostatic potential, z_{i} is the valence of species i, and ρ_{scale} is some characteristic number density. μ_{i}^{0} (x) is, for example, a hard-wall potential. With these definitions, Eq. (1) can be rewritten as a Poisson-drift-diffusion (Poisson-Nernst-Planck, PNP) system:

−J_{i} = D_{i} (x) A (x) [\frac{dρ_{i}}{dx} + \frac{1}{kT} ρ_{i} (x) (z_{i}e\frac{dϕ}{dx} + \frac{dμ_{i}^{0}}{dx} + \frac{dμ_{i}^{ex}}{dx} )]  (5)
\frac{dJ_{i}}{dx} = 0  (6)
−\frac{1}{A (x)} \frac{d}{dx} (ϵ (x) A (x) \frac{dϕ}{dx}) = e ∑_{i} z_{i}ρ_{i} (x)  (7)

where ϵ (x) is the local dielectric coefficient. In this paper we consider the special case

ϵ (x) = constant  (8)

so the Poisson equation we use is

−\frac{ϵ}{A (x)} \frac{d}{dx} (A (x) \frac{dϕ}{dx}) = e ∑_{i} z_{i}ρ_{i} (x) .  (9)

Position-dependent dielectric coefficients will be considered in future work.

4 Density functional theory

To calculate the excess chemical potentials μ_{i}^{ex} we use the density functional theory of Rosenfeld [6], [7]. We start by separating the grand potential into ideal (id), hard-sphere (HS), and electrostatic components (ES):

Ω ({ρ_{k} (x)}) = Ω_{id} ({ρ_{k} (x)}) + Ω_{HS} ({ρ_{k} (x)}) + Ω_{ES} ({ρ_{k} (x)})  (10)

with

Ω_{id} (\{ρ_{k} (x)\}) = kT ∑_{i} ∫ ρ_{i} (x) [ln (\frac{ρ_{i} (x)}{ρ_{scale}}) − 1 + \frac{1}{kT} [μ_{i}^{0} (x) − μ_{i} (x)]] dx.  (11)

The excess chemical potential μ_{i}^{ex} [Eq. (3)] is then the sum of the HS and ES components

μ_{i}^{ex} = μ_{i}^{HS} + μ_{i}^{ES}  (12)

where

μ_{i}^{HS} (x) = \frac{δΩ_{HS} (\{ρ_{k} (x′)\})}{δρ_{i} (x)}  (13)

and

μ_{i}^{ES} (x) = \frac{δΩ_{ES} (\{ρ_{k} (x′)\})}{δρ_{i} (x)}.  (14)

The DFT describes equilibrium systems and thus our use of it in a steady-state transport system requires the approximation of local equilibrium, an assumption validated by Frink, Thompson, and Salinger [17] who studied steady-state diffusion in hard-sphere fluids, as well as Eisenberg, Klosek, and Schuss [18] who found that, in the high-friction limit, velocity distributions of particles following Langevin dynamics were the sum of an equilibrium (Maxwellian) term and a term proportional to the flux.

4.1 Hard-sphere component

For the hard-sphere component, we use the “antisymmetrized” excess free energy density [8, Eq. (27)] because it is the best currently available for particles in confined geometries:

Φ_{HS} (\{n_{α} (x)\}) = −n_{0} ln (1 − n_{3}) + \frac{n_{1}n_{2} − n_{V 1}n_{V 2}}{1 − n_{3}} + \frac{n_{2}^{3}}{24π (1 − n_{3})^{2}} (1 − \frac{n_{V 2}n_{V 2}}{n_{2}^{2}} )^{3} (15)

with

Ω_{HS} ({c_{k} (x)}) = ∫ Φ_{HS} ({n_{α} (x′)}) dx′  (16)

where

n_{α} (x) = ∑_{i} ∫ ρ_{i} (x′) ω_{i}^{(α)} (x′ − x) dx′  (α = 0, 1, 2, 3, V 1, V 2)  (17)
ω_{i}^{(2)} (r) = δ (|r| − R_{i})  (18)
ω_{i}^{(3)} (r) = θ (|r| − R_{i})  (19)
ω_{i}^{(V 2)} (r) = \frac{r}{|r|}δ (|r| − R_{i})  (20)
4πR_{i}^{2}ω_{i}^{(0)} (r) = 4πR_{i}ω_{i}^{(1)} (r) = ω_{i}^{(2)} (r)  (21)
4πR_{i}ω_{i}^{(V 1)} (r) = ω_{i}^{(V 2)} (r)  (22)

and R_{i} is the radius of species i. δ is the Dirac delta function and θ is the unit step function, θ (x > 0) = 0 and θ (x ≤ 0) = 1. Then the hard-sphere component μ_{i}^{HS} of μ_{i}^{ex} is given by [6], [7]

μ_{i}^{HS} (x) = kT ∑_{α} ∫ \frac{∂Φ_{HS}}{∂n_{α}} (x′) ω_{i}^{(α)} (x − x′) dx′.  (23)

The integrals in Eqs. (17) and (23) simplify in the slab geometry. By analytically integrating over two components of the three-dimensional vector x and assuming that the functions ρ_{i} (x) and ∂Φ_{HS}/∂n_{α} (x) are constant over these two components, these integrals become [7]

n_{α} (x) = ∑_{i} ∫_{x−R_{i}}^{x+R_{i}} ρ_{i} (x′) W_{i}^{(α)} (x′ − x) dx′  (24)
μ_{i}^{HS} (x) = kT ∑_{α} ∫_{x−R_{i}}^{x+R_{i}} \frac{∂Φ_{HS}}{∂n_{α}} (x′) W_{i}^{(α)} (x − x′) dx′  (25)

where

W_{i}^{(2)} (r) = 2πR_{i}  (26)
W_{i}^{(3)} (r) = π (R_{i}^{2} − r^{2})  (27)
W_{i}^{(V 2)} (r) = 2πr (0, 0, 1)  (28)

4πR_{i}^{2}W_{i}^{(0)} (r) = 4πR_{i}W_{i}^{(1)} (r) = W_{i}^{(2)} (r)  (29)
4πR_{i}W_{i}^{(V 1)} (r) = W_{i}^{(V 2)} (r) .  (30)

4.2 Electrostatic component

For the electrostatic component μ_{i}^{ES} of μ_{i}^{ex} we use a generalization of Rosenfeld’s method [7] of expanding Ω_{ES} in a functional Taylor series in powers of

Δρ_{i} (x) = ρ_{i} (x) − ρ_{i}^{ref} (x)  (31)

where c_{i}^{ref} (x) is a given, nonconstant reference concentration profile. Then, up to second order,

Ω_{ES} ({ρ_{k} (x)}) ≈Ω_{ES} ({ρ_{k}^{ref} (x)}) + kT ∑_{i} ∫ μ_{i}^{ES} ({ρ_{k}^{ref} (x)}) Δρ_{i} (x) dx  (32)

− \frac{kT}{2} ∑_{i,j} ∫ ∫ c_{ij}^{(2),ES} (x, x′) Δρ_{i} (x) Δρ_{j} (x′) dxdx′

where c_{ij}^{(2),ES} is the electrostatic component of the second-order direct correlation function of the reference fluid [7]. Eq. (14) gives [7], [8]

μ_{i}^{ES} (x) = μ_{i}^{ES} ({ρ_{k}^{ref} (x)}) − z_{i}eϕ (x) − kT ∑_{j} ∫ c_{ij}^{(2),ES} (x, x′) Δρ_{j} (x′) dx′.  (33)

4.2.1 Reference fluid component

In Rosenfeld’s implementation [7] the reference fluid was the bulk (homogeneous) fluid that was in equilibrium with the inhomogeneous fluid. In the study of ion channels this will not suffice, however, as the regions of the system are very diverse. Channels are embedded in a membrane that separates two bulk fluids whose ionic concentrations are different, experimentally controlled, and can range from < 1 μM to 2 M. More importantly, inside the channel there is a very high concentration of charges from the amino acid residues of the channel protein. These attract equally high concentrations of counter charge with the total concentration of tens of molar. Thus, each section of the geometry has a different “intrinsic” reference fluid composition. This is the reason for choosing a non-uniform reference fluid.

To define the reference concentrations ρ_{k}^{ref} (x) for all species (permeating and confined) assume that some concentrations ρ_{k}^{∗} (x) and an estimate of the screening length at every point s^{∗} (x) are given. (These will later be determined by a self-consistency iteration.) From these given concentrationS we will construct a reference fluid that is charge-neutral and described by average, non-local densities.

To determine the ρ_{k}^{ref} (x), we start by defining

ρ̄_{k} (x) = Aρ_{k}^{∗} (x)  (z_{k} ≥ 0)  (34)

and

ρ̄_{k} (x) = ABρ_{k}^{∗} (x) (z_{k} < 0) (35)

with the constants A and B determined by requiring the {ρ̄_{k} (x)} to be charge-neutral and to have the same ionic strength at each location x as the {ρ_{k}^{∗} (x)}; that is, we require that

0 = ∑_{k} z_{k}ρ̄_{k} (x)  (36)

and

∑_{k} z_{k}^{2}ρ̄_{k} (x) = ∑_{k} z_{k}^{2}ρ_{k}^{∗} (x) .  (37)

These conditions give XXXX

A = \frac{∑_{k} z_{k}^{2}ρ_{k}^{∗} (x)}{∑_{z_{k}≥0} z_{k}^{2}ρ_{k}^{∗} (x) + B ∑_{z_{k}<0} z_{k}^{2}ρ_{k}^{∗} (x)}  (38)

and

B = \frac{∑_{z_{k}≥0} z_{k}ρ_{k}^{∗} (x)}{∑_{z_{k}<0} |z_{k}| ρ_{k}^{∗} (x)}.  (39)

The ρ_{k}^{ref} (x) are then defined as spatially averaged ρ̄_{k} (x):

ρ_{k}^{ref} (x) = ∫ ρ̄_{k} (x′) w (x′, x) dx′  (40)

where the weight function w is a normalized form of the volume weight ω_{i}^{(3)} (Eq. (19)):

w (x′, x) = \frac{θ (|x′ − x| − R_{filter} (x))}{\frac{4π}{3}R_{filter}^{3} (x)}.  (41)

The radius of the sphere over which we average is an approximation of the electrostatic length scale, namely the capacitance radius (that is, the ion radius plus the screening length):

R_{filter} (x) = s^{∗} (x) + \frac{∑_{k} ρ̄_{k} (x) R_{k}}{∑_{k} ρ̄_{k} (x)} .  (42)

The averaging produces a kind of non-local reference concentration, removing discontinuities that may have been in the {ρ̄_{k} (x)}. Futhermore, because we use the same filter and filter length for each species k, the resulting {ρ_{k}^{ref} (x)} are charge-neutral at every point since the {ρ̄_{k} (x)} are.

To calculate the ES excess chemical potential of the reference fluid, we apply the MSA for bulk fluids at every point [REFs]. This not only results in values for μ_{i}^{ES} ({ρ_{k}^{ref} (x)}), but also for the screening length, which we now redefine as

s (x) = \frac{1}{2Γ (x)}  (43)

where Γ (x) is the MSA parameter Γ that has been calculated at every point x. This then defines the capacitance length of each ion species:

λ_{k} (x) = R_{k} + λ (x) .  (44)

When compared against Monte Carlo simulations, the results obtained with reference fluid are very good after the system has been iterated to self-consistency (Sec. 4.4).

4.2.2 Non-reference fluid component

Calculating the last term of Eq. (33) requires the correlation function of the reference fluid c_{ij}^{(2),ES} (x, x′). We start by approximating the direct correlation function with the MSA [19]:

c_{ij}^{(2),ES} (x, x′) = −\frac{1}{kT} ψ_{ij} (x, x′)  (45)

if

|x − x′| > R_{ij} = R_{i} + R_{j}  (46)

where ψ_{ij} (x, x′) is the interaction potential of two point particles with charge z_{i}e and z_{j}e located at x and x′, respectively. In the case of uniform dielectric coefficient,

ψ_{ij} (x, x′) = \frac{z_{i}z_{j}e^{2}}{4πϵϵ_{0}} \frac{1}{|x − x′|}.  (47)

Since

kT ∑_{j} ∫ ψ_{ij} (x, x′) ρ_{j} (x′) dx′ = z_{i}eϕ (x) ,  (48)

we have

μ_{i}^{ES} (x) = μ_{i}^{ES} ({ρ_{k}^{ref} (x)})  (49)

− kT ∑_{j} ∫_{|x−x_{′}|≤R_{ij}} (c_{ij}^{(2),ES} (x, x′) + ψ_{ij} (x, x′)) Δρ_{j} (x′) dx′

where the core-overlap electrostatic correlations remain to be determined. A reasonable choice is the bulk MSA ES correlation function which can be expressed analytically even for mixtures [20], [21]. We choose not to use this formulation because the bulk MSA ES correlation function involves a uniform dielectric coefficient. While we use a uniform dielectric in this paper, to study ion permeation through ion channels more correctly, one must eventually use an inhomogeneous dielectric coeffient [14]. We will do this in future work and this will require an appropriate ES correlation function. The approximate ES correlation function we use in this paper is one that can be extended to the inhomogeneous dielectric coefficient case; we test it here for the uniform dielectric coefficient case, before applying it to more general situations. Specifically, we use the approximation given by Blum and Rosenfeld [21, Eq. (56)] that the ES correlation function c_{ij}^{(2),ES} is the negative of the interaction potential of two spherical shells Θ_{i} and Θ_{j} with each Θ_{k} having charge equal and opposite of the ion species k and radius λ_{k}, the capacitance radius of ion species k defined by Eq. (44).

With this ansatz, the short-range electrostatic correlations reduce to the interaction energy of the corresponding point charges z_{i}e and z_{j}e minus the interaction energy of particles of radius R_{i} and R_{j} with surface charge z_{i}e/4πR_{i}^{2} and z_{j}e/4πR_{j}^{2}, respectively [7], [8]:

kT c_{ij}^{(2),ES} (x, x′) = −\frac{z_{i}z_{j}e^{2}}{8πϵϵ_{0}} \frac{1}{|x − x′|} (|x − x′| \frac{λ_{i} + λ_{j}}{λ_{i}λ_{j}} − \frac{|x − x′|^{2}}{2λ_{i}λ_{j}} − \frac{(λ_{i} − λ_{j})^{2}}{2λ_{i}λ_{j}} ) (50)

In the slab geometry, Eq. (49) becomes

μ_{i}^{ES} (x) =μ_{i}^{ES} (\{ρ_{k}^{ref} (x)\}) − ∑_{j} \frac{z_{i}z_{j}e^{2}}{8ϵϵ_{0}} \frac{1}{λ_{i}λ_{j}}  (51)

× ∫_{x−R_{ij}}^{x+R_{ij}} c_{i} (x′) (\frac{1}{3}R_{ij}^{3} + λ_{ij}^{2}R_{ij} − λ_{ij}R_{ij}^{2} + λ_{ij}(x′ − x)^{2}

−\frac{1}{3} |x′ − x|^{3} − λ_{ij}^{2} |x′ − x|) dx′

where

λ_{ij} = λ_{i} + λ_{j}.  (52)

4.3 Iterating to self-consistency

I think this section is fabulous science and very good writing.

The choice of reference fluid is particularly important because it must be close enough to the final solution of the system that Δρ_{k} (x) is small so that the expansion of the ES grand potential (32) is valid. We achieve this by iteration on the choice of reference fluid. A guess of the solution serves as the initial set of {ρ_{k}^{∗} (x)} from which the reference concentrations are calculated as described in Sec. 4.2.1. The initial screening length s^{∗} (x) is set to zero everywhere. The system of equations is solved using this reference fluid and the updated screening and capacitance lengths defined in Eqs. (43) and (44), respectively. The concentrations from this solution then become the {ρ_{k}^{∗} (x)} for the next iteration cycle. After six iterations we observe no significant differences in the solutions.

4.4 Tests of the solution method

To ensure that our DFT method is valid, we compare our results for an equilibrium problem to Monte Carlo (MC) simulations. In the simulations, two compartments of different ionic compositions are equilibrated in the primitive model. In the left bath is a 24 M solution of half-charged oxygen ions that is confined to that compartment by a hard wall potential. In the right bath is a 0.1 M ionic solution. We conducted three separate tests with CaCl_{2}, NaCl, or KCl in the right compartment.

For these test simulations, the equations were solved as described below (Sec. 6) except that because the test XXX is an equilibrium problem we did not solve the Nernst-Planck equations (5), but rather

μ_{k} (x) = μ_{k}  (53)

for ion species where the total electrochemical potential μ_{k} was given for each ion. The results of the simulations (shown in Figs. 2-4) are quite good.

The particular ions and concentrations were chosen to approximate the calcium channel for which we later calculate current/voltage curves (Sec. 7). We vary the the cation species so that the accurancy we obtain bears directly on our ability to model ion selectivity.

5 Modeling the selectivity filter

Ion channels are proteins that can control permeation by placing some of their residues into their aqueous pore where they will XXXXX interact at close range with the permeating ions. It is these residues (many of which are charged) that confer the selectivity properties of the channel. COMMENT: It is important to use the same word for residues in both sentences. You can use amino acid if you wish but I prefer ‘residue’. For example, in the L-type calcium channel, there are four highlyconserved glutamate (E) residues (the “EEEE” locus) that produce the large physiological selectivity of calcium over sodium and potassium [REFs]. If these residues are mutated to aspartate (D), glutamate, lysine (K), and alanine (A) (the conserved “DEKA” locus of the voltage-gated sodium channel), then the channel becomes sodium selective under physiological conditions [REFs].

In equilibrium studies, these essential “structural charges” have been modeled as ions that are allowed to move freely inside the selectivity filter, but are confined there. They are not allowed to partition into the baths on each side of the channel [14], [MC refs]. We use the same approach here. Because these ions never reach the electrodes in the bath, they do not contribute to the measured current; each confined species j is in equilibrium:

\frac{dμ_{j}}{dx} = 0.  (54)

It is not a priori known what the value of the chemical potential μ_{j} is. What is known, however, is that the number of particles of each confined species j is fixed (N_{j}). Therefore, in conjunction with Eq. (54), for each confined species we solve the equation

∫_{x_{L}}^{x_{R}} ρ_{j} (x) dx = \frac{N_{j}}{V}  (55)

where the range of confinement is from x_{L} to x_{R} (the range of the selectivity filter) and V is the volume of the selectivity filter. COMMENT: I continue to think something should be said about this equation differing from the constant volume treatment of Hansen. Please explain to me if I am wrong, why I am wrong. If I am right, i.e., that this is very different in the two treatments, a remark should be made however politely. COMMENT: for next paragraph. I liked what I suggested better. I do not mind being over-ruled but please explain why. I give what I wrote again for this section below, separated by marker lines

Structural particles that are ions contribute to the electric field and therefore are also included in the Poisson equation (9) that calculates the electrostatic potential experienced by the permeating ions. However, because there is a small number of confined particles, if the electrostatic potential from Eq. (9) were used to calculate the electrochemical potential of the confined particles, then each confined particle would interact with a field that includes a contribution from the particle itself. This is a problem inherent to the mean-field approach, but one that is significant only when there are a small number OF species. COMMENT: I still do not understand the previous sentence except for the obvious that with an infinite number of particles NOT species, we approach the mean field. I think it essential to explain or reviewers will complain.

We alleviate this self-interaction by calculating the electrochemical potential of the confined particles with a different electrostatic potential for each confined species that applies only to that confined particle species. Specifically, in the calculation of μ_{j}, the electrostatic potential ϕ_{j} for confined species j is calculated from the modified Poisson equation

−\frac{ϵ}{A (x)} \frac{d}{dx} (A (x) \frac{dϕ_{j}}{dx} (x)) = e ∑_{k≠j} z_{k}ρ_{k} (x) + ez_{j} \frac{N_{j} − 1}{N_{j}}ρ_{j} (x)  (56)

where the concentration of confined species j is reduced by the factor 1/N_{j} to approximately eliminate the self-interaction.

######## PREVIOUS SUGGESTION FOLLOWS ###### comment: ions do exert force on themselves when there is a dielectric discon-

tinuity nearby. That is what the Appendix of EKS is all about. The key point is that ions should not exert force on themselves when they are in a medium of dielectric constant one, or in a medium with uniform dielectric constant. Now, if the dielectric fluctuates, and boundaries are created transiently, it is not clear what to do, but fortunately that is not our problem here (last sentence is a joke!!)

......

∫_{x_{L}}^{x_{R}} ρ_{j} (x) dx = \frac{N_{j}}{V}  (57)

where the range of confinement is from x_{L} to x_{R} (the range of the selectivity filter) and V is the volume of the selectivity filter.

This equation states the condition that is so different in our treatment compared to one that assumes a fixed structure and a rigid channel wall.

Because the structural particles are ions, they contribute to the electric field and therefore are also included in the Poisson equation (9) that calculates the electrostatic potential experienced by the permeating ions. However, because there is a finite number of confined particles, if the electrostatic potential from Eq. (9) were used to calculate the electrochemical potential of the confined particles, then each confined particle would interact with a field that was created by itself. The ion would create a force on itself, even if the domain were infinite without boundary conditions (REF Schuss, Nadler and Eisenberg, Appendix). This is a problem inherent to the mean-field approach, but one that is significant only when there are a finite number of ions of one species. The reader will want to know why this problem is only true when there a finite number. So do I! The approach we choose to alleviate this interaction is to allow the electrostatic force on each ion to differ, even if they have the same charge. We calculate the electrochemical potential of the confined particles with a different electrostatic potential for each confined species that applies only to that confined particle species. Specifically, in the calculation of μ_{j}, the electrostatic potential ϕ_{j} for confined species j is calculated from the modified Poisson equation

−\frac{ϵ}{A (x)} \frac{d}{dx} (A (x) \frac{dϕ_{j}}{dx} (x)) = e ∑_{i≠j} z_{i}c_{i} (x) + ez_{j} \frac{N_{j} − 1}{N_{j}}c_{j} (x)  (58)

where the concentration of confined species j is reduced by 1/N_{j} to approximately eliminate theself-interaction.

##################### END OF PREVIOUS SUGGESTION #########################

6 Numerical implementation

The inputs to the problem are the left and right far-bath concentrations of the permeating particles (ions and solvent), the applied electrostatic potential, and the number of structural (confined) ions, as well as the radii and diffusion coeffients of all particles and the dielectric coefficient. The outputs are the concentrations of all particles (permeating and structural) and the electrostatic potential, both as functions of location. These are found by simultaneously solving Eqs. (5), (6), (9), (25), (51), (54), (55), and (56). After these have been solved, fluxes are computed from an integrated Nernst-Planck equation (5):

J_{i} = \frac{ρ_{i} (L) exp (\frac{1}{kT} [μ_{i} (L) − ln (ρ_{i} (L))]) − ρ_{i} (R) exp (\frac{1}{kT} [μ_{i} (R) − ln (ρ_{i} (R))])}{∫ \frac{exp(\frac{1}{kT} [μ_{i}(x)−ln(ρ_{i}(x))])}{D_{i}(x)A(x)}dx}  (59)

where the integral is over the entire system and ρ_{i} (L/R) are the given left (L) and right (R) bath concentrations (similarly for the chemical potentials μ_{i} (L/R)).

6.1 Discretization

The grid we choose is nonuniform because the baths are necessarily much longer than the channel to which the structural charges are confined and, because the baths are of fixed concentration far from the channel, do not require the high density of grid points that the channel does where the functions are changing rapidly. In the channel we choose a grid with uniform spacing of 0.02 nm. (Numerical tests showed no significant differences when smaller spacings were used.) The spacing of the nonuniform grid outside the channel is chosen to be proportional to A (x), the area available for flux [4].

The differential equations (5), (6), (9), and (56) are rewritten as follows before they are discretized:

0 = \frac{d}{dx} [F_{i} (x) \frac{d}{dx} (u_{i} + \frac{z_{i}e}{kT} ϕ (x) + \frac{1}{kT} μ_{i}^{0} (x) + \frac{1}{kT} μ_{i}^{ex} (x))]  (60)
u_{i} (x) = ln [ρ_{i} (x)]  (61)

F_{i} (x) = \frac{D_{i} (x) A (x)}{min_{x} \{D_{i} (x) A (x)\}} exp [u_{i} (x)]  (62)
−\frac{ϵϵ_{0}kT}{N_{A}e^{2}d^{2}} \frac{1}{A (x)} \frac{d}{dx} (A (x) \frac{dϕ}{dx}) = ∑_{i} z_{i} exp [u_{i} (x)]  (63)

where all concentrations have been scaled by Avogadro’s number N_{A} and all lengths have been scaled by the system length d. (Eq. (56) may be rewritten in a similar form as Eq. (63)).

Both the recast Nernst-Planck equation (60) and the recast Poisson equation (63) have the same form for the derivatives:

\frac{d}{dx} (f (x) \frac{dy}{dx}) .  (64)

We discretize this derivative on the interior points of the grid {x_{0} = 0, x_{1}, x_{2}, …, x_{N−2}, x_{N−1} = 1} by

\frac{d}{dx} (f (x) \frac{dy}{dx})
≈\frac{2}{h_{m} + h_{m+1}} [f_{m+ \frac{1}{2}} \frac{dy}{dx} (x_{m+ \frac{1}{2}} ) − f_{m− \frac{1}{2}} \frac{dy}{dx} (x_{m− \frac{1}{2}} )]  (65)
≈\frac{2}{h_{m} + h_{m+1}} [\frac{f_{m} + f_{m+1}}{2} \frac{y_{m+1} − y_{m}}{h_{m+1}} − \frac{f_{m−1} + f_{m}}{2} \frac{y_{m} − y_{m−1}}{h_{m}} ]  (66)
= (f_{m} + f_{m+1}) (y_{m+1} − y_{m}) β_{1,m} − (f_{m−1} + f_{m}) (y_{m} − y_{m−1}) β_{−1,m}  (67)

where

f_{m} = f (x_{m})  (68)
y_{m} = y (x_{m})  (69)
β_{m,−1} = \frac{1}{h_{m} (h_{m+1} + h_{m})} β_{m,1} = \frac{1}{h_{m+1} (h_{m+1} + h_{m})} (70)
h_{m} = x_{m} − x_{m−1}.  (71)

This is finally rewritten as

\frac{d}{dx} (f (x) \frac{dy}{dx})
≈β_{−1,m} (f_{m−1} + f_{m}) y_{m−1} − [β_{−1,m} (f_{m−1} + f_{m}) + β_{1,m} (f_{m} + f_{m+1})] y_{m}  (72)
+ β_{1,m} (f_{m} + f_{m+1}) y_{m+1}.

All functions are specified on the boundary nodes x_{0} and x_{N−1}, either because the function has been specified for all grid points (the diffusion coefficients, for example) or because of boundary conditions for the problem (the electrostatic potential, for example). Thus these derivatives are defined on the interior grid points, x_{1}, …, x_{N−1}.

The integrals in Eqs. (24), (25), and (51) have a similar structure and thus the integration scheme we choose was the same for all three. These integrals are over the range x−r to x+r with different r’s for each integral. x is always a grid point, but r is a radius and because of the nonuniform grid and different ion sizes, it not generally true that x − r or x + r are grid points. Thus it is necessary to be able to evaluate these integrals at points between two grid nodes. We do this by dividing the integral into a sum of integrals, with each new integral ranging over two consecutive grid points or the last grid point to the endpoint of the original integral. Since in Eqs. (24), (25), and (51) one of the functions can be written as a polynomial in the integration variable x′, each new integral in the sum can be evaluated analytically by assuming that the other function in the integral (ρ_{i} (x′) or ∂Φ_{HS}/∂n_{α} (x′)) is a linear function in x′ between consecutive grid nodes.

The remaining equations to be discretized are Eqs. (54) and (55) for the confined particles. Let x_{M_{L}} and x_{M_{R}} be the left and right grid nodes, respectively, on and between which the structural particles are confined. Then Eq. (54) is discretized as

0 = μ_{j} (x_{m}) − μ_{j} (x_{m+1})  (73)

for M_{L} ≤ m ≤ M_{R} − 1. For the last grid node of confinement x_{M_{R}}, we apply Eq. (55). Because this integral is always evaluated only in the region where the grid spacing is uniform and small, we evaluate this with the one-half rule. COMMENT: readers may wonder, as I do, what is the "one-half rule"??

6.2 Solving the discretized system

After the system of equations is discretized, the values of u_{i} (x_{m}) = ln [ρ_{i} (x_{m})] for all permeant species i and ϕ (x_{m}) remain to be determined at all interior grid points (m = 1, …, N−2), as well as u_{j} (x_{m}) = ln [ρ_{i} (x_{m})] and ϕ_{j} (x_{m}) for all confined species j at all grid points to which they are confined (m = M_{L}, …, M_{R}). We solved the system of discretized equations using Newton’s method [23] because of its square-convergence properties, because it is possible to analytically evaluate the Jacobian (derivative) matrix of the discretized equations, and, most importantly, because iterative methods that did not update all variables at every iteration did not converge, COMMENT: comma was inserted or required a very large number of iterations. Using Newton’s method with a poor initial guess, 30 iterations are typical with a total calculation time of one to two minutes, depending on the number ion species. Many fewer iterations (4-10) are necessary when a better inital guess is available, for example when calculating a current/voltage curve. Then a good intial guess is a previous solution of the system with a different applied voltage (up to approximately 0.1 V difference).

Because of the specific functional dependence of the equations, it is possible to write the Jacobian matrix as a band matrix. Furthermore, without the discretized version of Eq. (55) the band becomes significantly more narrow. Because Eq. (55) is a solution condition for Eq. (54) and is only used once for each confined species, it was more efficient to solve the Jacobian matrix equation for each Newton iteration in two steps, first solving the narrow band matrix without Eq. (55) using an LU decomposition for band matrices [23] and finally solving the complete matrix equation with the Sherman-Morrison-Woodbury method [23]. The Newton iteration was stopped when the absolute value of the largest difference between the same variable from two consecutive iterations was less than 10^{−8}.

7 Example: a calcium channel

To give an example of the outputs of our algorithm, we model a calcium-selective channel. In this channel, the structural particles are four fully-charged glutamate residues as those found in the L-type calcium channel [4], [14]. We model each glutamate residue by the carboxyl group at its end. Each of these carboxyls we model as two unbonded, half-charged oxygen ions (diameter 0.28 nm) that are confined with a hard-wall potential to the cylindrical selectivity filter (channel) of radius 0.35 nm and length 1 nm. Throughout the system we assume that the dielectric coefficient is uniform with ϵ = 78.4, the value for pure water. Previous studies have indicated that inside the channel an effective dielectric coefficient of approximately 10 is needed to have the channel select calcium in the micromolar range as the real protein do [14]; the choice of a uniform dielectric coefficient will cause our simulated channel to be less selective than the real channel. Furthermore, these studies showed that including water as uncharged, hard spheres gives significantly different results than using the primitive model and therefore in the current simulations we also include water in this way. The values of the other input parameters are listed in the legend of Fig. 5.

For this channel we simulate the classical selectivity experiment of the L-type calcium channel where the left and right baths contain 100 mM NaCl and CaCl_{2} is added to the right bath. In Fig. 5 are shown the concentration profiles of Ca^{2+}, Na^{+}, Cl^{−}, and O^{−1/2} for different CaCl_{2} concentrations in the right bath. As CaCl_{2} is added, Na^{+} is displaced from the channel, with more than half of the Na^{+} displaced by ?? of CaCl_{2}, showing that the channel is calcium selective. This can also be seeing in the current/voltage curves for this experiment (shown in Fig. 6); at ?? of CaCl_{2} the conductance at 0 mV is half of the two limiting conductances. For the natural L-type calcium channel, this occurs at 1 μM, indicating that while our model gives significant calcium selectivity, other aspects are needed to fully model the natural channel.

Acknowledgments

This work was supported by grants from DARPA (R.S.E. and W.N.) and NIH (T32NS07044 to D.G.). We greatfully acknowledge Dezso Boda for providing the Monte Carlo simulations used in this paper.

References

[1] B. Hille, Ion Channels of Excitable Membranes, 3rd ed. (Sinauer Associates, Inc., Sunderland, MA, 2001).

[2] D. P. Chen, J. Lear, and B. Eisenberg, Biophys. J. 72, 97 (1997).

[3] D. P. Chen, et al., Biophys. J. 76, 1346 (1999).

[4] W. Nonner and B. Eisenberg, Biophys. J. 75, 1287 (1998).

[5] C. Caccamo, J.-P. Hansen, and G. Stell, New Approaches to Problems in Liquid State Theory: Inhomogeneities and Phase Separation in Simple, Complex, and Quantum Fluids (Kluwer Academic Publishers, Dordrecht, The Netherlands, 1999).

[6] Y. Rosenfeld, Phys. Rev. Lett. 63, 980 (1989).

[7] Y. Rosenfeld, J. Chem. Phys. 98, 8126 (1993).

[8] Y. Rosenfeld, M. Schmidt, H. Löwen, and P. Tarazona, Phys. Rev. E 55, 4245 (1997).

[9] T. Biben, J.-P. Hansen, and Y. Rosenfeld, Phys. Rev. E 57, R3727 (1998).

[10] B. Sakmann and E. Neher, eds., Single-Channel Recording (Plenum Press, New York, 1983).

[11] O. P. Hamill, A. Marty, E. Neher, B. Sakmann, and F. J. Sigworth, Pflügers Arch. 391, 85 (1981).

[12] J. D. Romano and R. H. Price, Am. J. Phys. 64, 1150 (1996).

[13] D. Goulding, S. Melchionna, and J.-P. Hansen, Phys. Chem. Chem. Phys. 3, 1644 (2001).

[14] W. Nonner, D. Gillespie, D. Henderson, and B. Eisenberg, J. Phys. Chem. B 105, 6427 (2001).

[15] B. Corry, S. Keyucak, and S.-H. Chung, Biophys. J. 78, 2364 (2000).

[16] Z. Schuss, B. Nadler, and B. Eisenberg, Phys. Rev. E 64, 036116 (2001).

[17] L. J. D. Frink, A. Thompson, and A. G. Salinger, J. Chem. Phys. 112, 7564 (2000).

[18] R. S. Eisenberg, M. M. Klosek, and Z. Schuss, J. Chem. Phys. 102, 1767 (1995).

[19] E. Waisman and J. L. Lebowitz, J. Chem. Phys. 56, 3086 (1972).

[20] K. Hiroike, Mol. Phys. 33, 1195 (1977).

[21] L. Blum and Y. Rosenfeld, J. Stat. Phys. 63, 1177 (1991).

[22] D. Gillespie, Ph.D. thesis, Rush University, Chicago (1999).

[23] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, 2nd ed. (Cambridge University Press, Cambridge, UK, 1992).