Coupling Poisson-Nernst-Planck and Density

Functional Theory to Calculate Ion Flux

Wolfgang Nonner  Dirk Gillespie  Robert S. Eisenberg

May 26, 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, often charged amino acids. Enough data on 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. XXXXX 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 spherical 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 charged spheres confined to the channel. In this one dimensional model (unlike some others we have studied PNP references here), 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 [2], [3], [4], [5], [6]. 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

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 experiments [7], [8]. comment: what a great sentence (the next one; last one is good too)! 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 includes the effects of the finite volume of ions and functional groups of the protein and 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 the baths: 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. 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 [9].

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 hard wall at the pore/protein interface that maintained a fixed volume for the pore, large radial packing effects would result in the pore [10]. However, real proteins are not rigid and therefore we take the opposite view of a hard-wall interface: In our model, the interface separates regions of different pressure. 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 locatoin, but rather by an average cross-section area or pore volume that may be constant or variable depending on the compressability of the protein [11]. Here we consider the simple case that the protein maintains a constant average pore volume. comment: the reader will be confused by the last sentence. If we are maintaining constant volume, isn’t this equivalent to the Hansen treatment. It is ESSENTIAL to explain the difference explicitly. By the way the previous paragraph is very very good and goes much further to explain this subtle idea than in any of our other papers. Note I moved the next paragraph.

It is important to note that this representation of the pore/protein interface is appropriate even necessary (in our opinion) even in three dimensions. Hard walls are unnatural when describing protein channels and are likely to introduce behavior not found in systems with flexible walls. It is not introduced to make the problem one dimensional.

Using this description of the interface, we can 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): Comment: the logical connection between the parts of this sentence are not clear. Why does the constant pressure description make this description, e.g., slab, possible?? This is a presentation not scientific issue, of course.

1. From the hard-sphere perspective, the protein makes the pore effectively wider than its physical radius; that is, seen as only a hard-sphere liquid, the combination of the pore and protein is approximately one continuous liquid and because it is relatively large, it may be approximated in the slab geometry. Comment: please explain. I do not understand and perhaps the reader will have the same problem.

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 many cases (such as the example calcium channel we will consider later), there are such a large concentrations of ions inside the channel that the Debye length less than half the physical radius of the pore. Brownian dynamics studies have shown that under these conditions slab geometry can be used in mean-field theories like the one we consider here [12].

3 Modeling ion flux

We describe steady-state flux by the constituitive law (Ohm’s or Fick’s law for electrolytes)

−J_{i} = \frac{1}{kT} D_{i} (x) A (x) c_{i} (x) \frac{dμ_{i}}{dx} (x)  (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} and the gradient of the electrochemical potential μ_{i} (x). A (x) is the area of the spherical shells described in the previous section.

The electrochemical potential consists of the ideal component μ_{i}^{id} (x), the excess component μ_{i}^{ex} (x), 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 [c_{i} (x)]  (4)

where ϕ (x) is the local electrostatic potential and z_{i} is the valence of species i. μ_{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{dc_{i}}{dx} (x) + \frac{z_{i}e}{kT} c_{i} (x) \frac{dϕ}{dx} (x) + \frac{1}{kT} c_{i} (x) \frac{dμ_{i}^{0}}{dx} (x) + \frac{1}{kT} c_{i} (x) \frac{dμ_{i}^{ex}}{dx} (x))

(5)

\frac{dJ_{i}}{dx} = 0  (6)
−\frac{1}{A (x)} \frac{d}{dx} (ϵ (x) A (x) \frac{dϕ}{dx} (x)) = e ∑_{i} z_{i}c_{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} (x)) = e ∑_{i} z_{i}c_{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 [3], [4]. We start by separating the grand potential into ideal (id), hard-sphere (HS), and electrostatic components (ES):

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

with

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

where λ_{i} is the deBroglie wavelength. The deBroglie wavelength is [physical phrase or clause] and is given by [insert equation].

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} (\{c_{k} (x′)\})}{δc_{i} (x)}  (13)

and

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

We note that DFT was designed for equilibrium systems and thus our use of it in a steady-state transport system requires the assumption of local equilibrium, an assumption validated by Eisenberg, Klosek and Schuss, REF who studied systems that might not be in local equilibrium, and by Frink, Thompson, and Salinger [14] who studied the steady-state diffusion of colored hard spheres.

4.1 Hard-sphere component

For the hard-sphere component, we use the “antisymmetrized” excess free energy density [5, Eq. (27)] Comment: The reader will want to know why you chose this.

Φ_{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} ∫ c_{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 [3], [4]

μ_{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 c_{i} (x) and ∂Φ_{HS}/∂n_{α} (x) are constant over these two components, these integrals become [4]

n_{α} (x) = ∑_{i} ∫_{x−R_{i}}^{x+R_{i}} c_{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 [4] of expanding Ω_{ES} in a functional Taylor series in powers of

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

where c_{i}^{ref} (x) is a given, nonconstant reference concentration profile. (For simplicity, we choose this reference fluid to be charge-neutral everywhere.) Then, up to second order,

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

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

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

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

4.2.1 Non-reference fluid component

To approximate the last term of Eq. (33) we use the method of Biben, Hansen, and Rosenfeld [5] which applies the Mean Spherical Approximation (MSA) in the strong coupling limit. Specifically, when two ions are far apart in this limit, they interact only through the electric mean field and other correlations become insignificant.The MSA approximates all correlations outside of the contact radius

R_{ij} = R_{i} + R_{j}  (34)

by this asymptotic limit [Waisman&Lebowitz]. In our case this gives

kT ∑_{j} ∫_{|x−x_{′}|>R_{ij}} c_{ij}^{(2),ES} (x, x′) c_{j} (x′) dx′ ≈ −z_{i}eϕ (x)  (35)

or [5]

μ_{i}^{ES} (x) = μ_{i}^{ES} ({c_{k}^{ref} (x)}) − kT ∑_{j} ∫_{|x−x_{′}|≤R_{ij}} c_{ij}^{(2),ES} (x, x′) Δc_{j} (x′) dx′.  (36)

Lastly, 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 [4], [5]:

kT c_{ij}^{(2),ES} (x, x′)

= \frac{z_{i}z_{j}e^{2}}{4πϵ (x) ϵ_{0}} \frac{1}{|x − x′|} − \frac{z_{i}z_{j}e^{2}}{8πϵ (x) ϵ_{0}} \frac{1}{|x − x′|} (|x − x′| \frac{R_{i} + R_{j}}{R_{i}R_{j}} − \frac{|x − x′|^{2}}{2R_{i}R_{j}} − \frac{(R_{i} − R_{j})^{2}}{2R_{i}R_{j}} )

(37)

= \frac{z_{i}z_{j}e^{2}}{16πϵ (x) ϵ_{0}} \frac{1}{R_{i}R_{j}} \frac{1}{|x − x′|} [|x − x′| − (R_{i} + R_{j})]^{2}  (38)

where we have used the simplifying assumption that the dielectric coefficient is approximately constant with a value of ϵ (x). This is exact in the case of a constant dielectric coefficient ϵ. Comment: the previous sentence needs to be explained. Why is it exact then? In the slab geometry, Eq. (36) becomes

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

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

4.2.2 Reference fluid component

The only remaining component of the excess chemical potential is the electrostatic excess chemical potential of the reference fluid, μ_{i}^{ES} ({c_{k}^{ref} (x)}) in Eq. (39). In Rosenfeld’s implementation [4] 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 geometry is very different. Channels are embedded in a membrane that separates two bulk fluids whose concentrations are different, experimentally controlled, and can range from 1 nM to 2 M. Furthermore, inside the channel there is a 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, the channel geometry has three different “intrinsic” reference fluid compositions. This is the reason for choosing a non-constant reference fluid.

Comment: The previous paragraph is a brilliant piece of science. Congratulations!!!!!!!!!

The reference fluid we choose is divided into three geometrically distinct regions: the left bath, channel, and right bath. The left and right bath concentrations of all ions are experimentally controlled and thus we choose these as our reference fluids in their respective geometrical regions. Choosing the reference fluid in the channel region is more difficult because the ion concentrations inside the channel are not known until after the analysis which depends on a reference state that itself depends on the ion concentrations. Ion concentrations are an output of the calculations. We do, however, have an understanding of the approximate ion concentrations inside the channel under some experimental conditions, and these are what we choose to use as the channel reference fluid. Having chosen the piecewise-constant compositions, the cation concentrations are Gaussfiltered Comment: Reference is needed here to some book on Gaussian filters. Also say why you chose this filter characteristic, presumably to change time dependent functions as little as possible to produce continuous concentration profiles. A charge-neutral reference fluid is produced by adjusting the anion concentrations.

Scientific Comment: for future work AFTER this paper ONLY. If the dependence of reference state on concentration is know by formula or simple computer program, one might be able to do the calculation selfconsistently. That is assume the above reference state, calculate the concentrations, calculate a new reference state, etc.

The resulting charge-neutral reference fluid is constant in the three geometric regions and smoothly connected by the Gauss-filtering described above. Each constant region is then a homogeneous fluid which we describe with the MSA [11]. The connecting regions are not homogeneous fluids and to describe them correctly would require a (unavailable) theory of inhomogeneous fluids. However, because these regions are connected to two homogeneous fluids, the correct electrostatic excess chemical potential is some interpolation of the excess chemical potentials of the two bulk fluids. To approximate the electrostatic excess chemical potential in the connecting regions, we treat each point is a bulk liquid and again use the MSA to calculate the chemical potentials, including its component the excess electrostatic excess chemical potential.

5 Modeling the selectivity filter

Ion channels are proteins and some of the amino acid residues can enter the permeation pathway and interact with the permeating particles. It is these residues (many of which are charged) that confer the selectivity properties of the channel. For example, in the L-type calcium channel, there are four highly-conserved glutamate (E) residues (the “EEEE” locus) that produce the large selectivity of sodium over calcium [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 [REFs].

In equilibrium studies, these “structural charges” have been modeled as spherical charges (ions) that are allowed to move freely inside the selectivity filter, but are not allowed to partition into the baths on each side of the channel [11], [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.  (40)

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

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

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)  (42)

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

Comment: this is a brilliant solution to the problem. Now let’s see Zeev prove it is right!!!

6 Numerical implementation

The inputs to the problem are the left and right bath concentrations of the permeating ions, 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), (39), (40), (41), and (42).

6.1 Discretization

The grid we chose was nonuniform because the baths are necessarily much longer than the channel to which the structural charges are confined and, because XXX 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 difference when smaller spacings were used.) The spacing of the nonuniform grid outside the channel was chosen to be proportional to A (x), the area available for flux [16].

The differential equations XXX (5), (6), (9), and (42), 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))]  (43)
u_{i} (x) = ln [c_{i} (x)]  (44)

F_{i} (x) = \frac{D_{i} (x) A (x)}{min_{x} \{D_{i} (x) A (x)\}} exp [u_{i} (x)]  (45)
−\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)]  (46)

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

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

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

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}} )]  (48)
≈ \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}} ]  (49)
= (f_{m} + f_{m+1}) (y_{m+1} − y_{m}) β_{1,m} − (f_{m−1} + f_{m}) (y_{m} − y_{m−1}) β_{−1,m}  (50)

where

f_{m} = f (x_{m})  (51)
y_{m} = y (x_{m})  (52)
β_{m,−1} = \frac{1}{h_{m} (h_{m+1} + h_{m})} β_{m,1} = \frac{1}{h_{m+1} (h_{m+1} + h_{m})} (53)
h_{m} = x_{m} − x_{m−1}.  (54)

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}  (55)
+ β_{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 (39) XXXX 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 (39) 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 (c_{i} (x′) or ∂Φ_{HS}/∂n_{α} (x′)) is a linear function in x′ between consecutive grid nodes.

The remaining equations to be discretized are Eqs. (40) and (41) 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. (40) is discretized as

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

for M_{L} ≤ m ≤ M_{R} − 1. For the last grid node of confinement x_{M_{R}}, we apply Eq. (41). Because this integral is always evaluated only in the region where the grid spacing is uniform and small, we evaluate this with a simple rectangular integration scheme.

6.2 Solving the discretized system

After XXX the system of equations is discretized, the values of u_{i} (x_{m}) = ln [c_{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 [c_{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 [13] 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 only converged (over the necessary domain of parameters) if they update all variables at every iteration.

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. (41) the band becomes significantly more narrow. Because Eq. (41) is a solution condition for Eq. (40) 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. (41) using an LU decomposition for band matrices [13] and finally solving the complete matrix equation with the Sherman-Morrison-Woodbury method [13]. 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}.

Comment: Please give an estimate of the time the calculation takes.

7 Example: a calcium channel

In order to give an example of the outputs of our algorithm and to explain some details of the implementation, we give a simple model of a calcium-selective channel. In such a channel, the structural particles are four fully-charge glutamate residues. 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.

References

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

[2] 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).

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

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

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

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

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

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

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

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

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

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

[13] 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).

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

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

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