%% This document created by Scientific Word (R) Version 3.5

\documentclass[12pt]{article}%
\usepackage{amsmath}
\usepackage{graphicx}%
\usepackage{amsfonts}%
\usepackage{amssymb}
%TCIDATA{OutputFilter=latex2.dll}
%TCIDATA{CSTFile=LaTeX article (bright).cst}
%TCIDATA{Created=Thursday, March 29, 2001 22:48:04}
%TCIDATA{LastRevised=Friday, May 24, 2002 13:52:56}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{<META NAME="DocumentShell" CONTENT="General\SW\Blank - Standard LaTeX Article">}
\newtheorem{theorem}{Theorem}
\newtheorem{acknowledgement}[theorem]{Acknowledgement}
\newtheorem{algorithm}[theorem]{Algorithm}
\newtheorem{axiom}[theorem]{Axiom}
\newtheorem{case}[theorem]{Case}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{conclusion}[theorem]{Conclusion}
\newtheorem{condition}[theorem]{Condition}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{criterion}[theorem]{Criterion}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{exercise}[theorem]{Exercise}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{notation}[theorem]{Notation}
\newtheorem{problem}[theorem]{Problem}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{solution}[theorem]{Solution}
\newtheorem{summary}[theorem]{Summary}
\newenvironment{proof}[1][Proof]{\textbf{#1.} }{\ \rule{0.5em}{0.5em}}
\textwidth 6.5in
\textheight 8.5 in
\oddsidemargin 0in
\topmargin -0.35 in
\renewcommand{\baselinestretch}{1.6}


\begin{document}

\title{Coupling Poisson-Nernst-Planck and Density Functional Theory to Calculate Ion Flux}
\author{Wolfgang Nonner
\and Dirk Gillespie
\and Robert S. Eisenberg}
\maketitle

\section{Introduction}

The movement of ions through open ion channels is one of the many interesting
physical chemistry problems provided by living cells. In this paper we
describe an 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 virtually
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 \cite{Hille}.

Enough data on functional groups (often charged amino acid residues) 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. Towards these goals, in this paper we describe
an engineering approach to selective ion conduction through ion channels.
Specifically, we develop a one-dimensional model of ion transport 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 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 \cite{Caccamo},
\cite{Rosenfeld89}, \cite{Rosenfeld93}, \cite{Rosenfeld97}, \cite{Rosenfeld98}%
. 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.

\section{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
\cite{SakmannNeher}, \cite{originalSN}. 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 limited 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. In this section we derive an approximation
that is fine grained in the flux-limiting region (the channel) and course
grained in the baths. This 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\left(  x\right)  $ 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
\cite{Romano}.

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 large radial
packing effects would result in the pore \cite{Hansen}. However, real proteins
are not 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, but rather by an average cross-section area or pore volume
that may be constant or variable depending on the compressability of the
protein \cite{NGHE}. Here we consider the simple case that the protein
maintains a constant average pore volume.

With this interface, it is possible to 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):

\begin{enumerate}
\item 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.

\item 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 will
become valid. 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 out 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.

\item 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 this scale difference is sufficient to apply
the slab geometry\ in mean-field theories like the one we consider here
\cite{Chung}.
\end{enumerate}

It is important to note that this representation of the pore/protein interface
is not meant as a means to a one-dimensional approximation; even in a
three-dimensional calculation a model of the protein like this would be needed
to correctly model the channel (hard walls are unnatural when discussing channels).

\section{Modeling ion flux}

We describe steady-state flux by the constituitive law%
\begin{equation}
-J_{i}=\frac{1}{kT}D_{i}\left(  x\right)  A\left(  x\right)  c_{i}\left(
x\right)  \frac{d\mu_{i}}{dx}\left(  x\right)  \label{NP_eqn_alt}%
\end{equation}%
\begin{equation}
\frac{dJ_{i}}{dx}=0;
\end{equation}
that is, for particle species $i$, the (constant) flux $J_{i}$ is proportional
to the particle density $\rho_{i}$ and the gradient of the electrochemical
potential $\mu_{i}\left(  x\right)  $. $A\left(  x\right)  $ is the area of
the spherical shells described in the previous section.

The electrochemical potential consists of the ideal component $\mu
_{i}^{\text{id}}\left(  x\right)  $, the excess component $\mu_{i}^{\text{ex}%
}\left(  x\right)  $, and the concentration-independent component $\mu_{i}%
^{0}\left(  x\right)  $:%
\begin{equation}
\mu_{i}\left(  x\right)  =\mu_{i}^{0}\left(  x\right)  +\mu_{i}^{\text{id}%
}\left(  x\right)  +\mu_{i}^{\text{ex}}\left(  x\right)  \label{chem_potl}%
\end{equation}
with%
\begin{equation}
\mu_{i}^{\text{id}}\left(  x\right)  =z_{i}e\phi\left(  x\right)
+kT\ln\left[  c_{i}\left(  x\right)  \right]  \label{ideal_def}%
\end{equation}
where $\phi\left(  x\right)  $ is the local electrostatic potential and
$z_{i}$ is the valence of species $i$. $\mu_{i}^{0}\left(  x\right)  $ is, for
example, a hard-wall potential. With these definitions, Eq. $\left(
\ref{NP_eqn_alt}\right)  $ can be rewritten as a Poisson-drift-diffusion
(Poisson-Nernst-Planck) system:
\begin{equation}
-J_{i}=D_{i}\left(  x\right)  A\left(  x\right)  \left(  \frac{dc_{i}}%
{dx}\left(  x\right)  +\frac{z_{i}e}{kT}c_{i}\left(  x\right)  \frac{d\phi
}{dx}\left(  x\right)  +\frac{1}{kT}c_{i}\left(  x\right)  \frac{d\mu_{i}^{0}%
}{dx}\left(  x\right)  +\frac{1}{kT}c_{i}\left(  x\right)  \frac{d\mu
_{i}^{\text{ex}}}{dx}\left(  x\right)  \right)  \label{NP_eqn}%
\end{equation}%
\begin{equation}
\frac{dJ_{i}}{dx}=0\label{continuity_eqn}%
\end{equation}%
\begin{equation}
-\frac{1}{A\left(  x\right)  }\frac{d}{dx}\left(  \epsilon\left(  x\right)
A\left(  x\right)  \frac{d\phi}{dx}\left(  x\right)  \right)  =e\sum_{i}%
z_{i}c_{i}\left(  x\right)  \label{gen_Poisson_eqn}%
\end{equation}
where $\epsilon\left(  x\right)  $ is the local dielectric coefficient. In
this paper we consider the special case%
\begin{equation}
\epsilon\left(  x\right)  =\text{constant}%
\end{equation}
so the Poisson equation we use is%
\begin{equation}
-\frac{\epsilon}{A\left(  x\right)  }\frac{d}{dx}\left(  A\left(  x\right)
\frac{d\phi}{dx}\left(  x\right)  \right)  =e\sum_{i}z_{i}c_{i}\left(
x\right)  .\label{Poisson_eqn}%
\end{equation}
Position-dependent dielectric coefficients will be considered in future work.

\section{Density functional theory}

To calculate the excess chemical potentials $\mu_{i}^{\text{ex}}$ we use the
density functional theory of Rosenfeld \cite{Rosenfeld89}, \cite{Rosenfeld93}.
We start by separating the grand potential into ideal (id), hard-sphere (HS),
and electrostatic components (ES):
\begin{equation}
\Omega\left(  \left\{  c_{k}\left(  \mathbf{x}\right)  \right\}  \right)
=\Omega_{\text{id}}\left(  \left\{  c_{k}\left(  \mathbf{x}\right)  \right\}
\right)  +\Omega_{\text{HS}}\left(  \left\{  c_{k}\left(  \mathbf{x}\right)
\right\}  \right)  +\Omega_{\text{ES}}\left(  \left\{  c_{k}\left(
\mathbf{x}\right)  \right\}  \right)  \label{grand_potl}%
\end{equation}
with%
\begin{equation}
\Omega_{\text{id}}\left(  \left\{  c_{k}\left(  \mathbf{x}\right)  \right\}
\right)  =kT\sum_{i}\int c_{i}\left(  \mathbf{x}\right)  \left\{  \ln\left[
\lambda_{i}^{3}c_{i}\left(  \mathbf{x}\right)  \right]  -1+\frac{1}{kT}\left[
\mu_{i}^{0}\left(  \mathbf{x}\right)  -\mu_{i}\left(  \mathbf{x}\right)
\right]  \right\}  d\mathbf{x}%
\end{equation}
where $\lambda_{i}$ is the deBroglie wavelength. The excess chemical potential
$\mu_{i}^{\text{ex}}$ [Eq. $\left(  \ref{chem_potl}\right)  $] is then the sum
of the HS and ES components%
\begin{equation}
\mu_{i}^{\text{ex}}=\mu_{i}^{\text{HS}}+\mu_{i}^{\text{ES}}%
\end{equation}
where%
\begin{equation}
\mu_{i}^{\text{HS}}\left(  \mathbf{x}\right)  =\frac{\delta\Omega_{\text{HS}%
}\left(  \left\{  c_{k}\left(  \mathbf{x}^{\prime}\right)  \right\}  \right)
}{\delta c_{i}\left(  \mathbf{x}\right)  }%
\end{equation}
and%
\begin{equation}
\mu_{i}^{\text{ES}}\left(  \mathbf{x}\right)  =\frac{\delta\Omega_{\text{ES}%
}\left(  \left\{  c_{k}\left(  \mathbf{x}^{\prime}\right)  \right\}  \right)
}{\delta c_{i}\left(  \mathbf{x}\right)  }. \label{mu_es_from_gp}%
\end{equation}

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 Frink, Thompson, and Salinger
\cite{frink} who studied the steady-state diffusion of colored hard spheres.

\subsection{Hard-sphere component}

For the hard-sphere component, we use the ``antisymmetrized'' excess free
energy density \cite[Eq. (27)]{Rosenfeld97}
\begin{equation}
\Phi_{\text{HS}}\left(  \left\{  n_{\alpha}\left(  \mathbf{x}\right)
\right\}  \right)  =-n_{0}\ln\left(  1-n_{3}\right)  +\frac{n_{1}n_{2}%
-n_{V1}n_{V2}}{1-n_{3}}+\frac{n_{2}^{3}}{24\pi\left(  1-n_{3}\right)  ^{2}%
}\left(  1-\frac{n_{V2}n_{V2}}{n_{2}^{2}}\right)  ^{3}%
\end{equation}
with%
\begin{equation}
\Omega_{\text{HS}}\left(  \left\{  c_{k}\left(  \mathbf{x}\right)  \right\}
\right)  =\int\Phi_{\text{HS}}\left(  \left\{  n_{\alpha}\left(
\mathbf{x}^{\prime}\right)  \right\}  \right)  d\mathbf{x}^{\prime}%
\end{equation}
where%
\begin{equation}
n_{\alpha}\left(  \mathbf{x}\right)  =\sum_{i}\int c_{i}\left(  \mathbf{x}%
^{\prime}\right)  \omega_{i}^{\left(  \alpha\right)  }\left(  \mathbf{x}%
^{\prime}-\mathbf{x}\right)  d\mathbf{x}^{\prime}\qquad\left(  \alpha
=0,1,2,3,V1,V2\right)  \label{n_def}%
\end{equation}%
\begin{equation}
\omega_{i}^{\left(  2\right)  }\left(  \mathbf{r}\right)  =\delta\left(
\left|  \mathbf{r}\right|  -R_{i}\right)
\end{equation}%
\begin{equation}
\omega_{i}^{\left(  3\right)  }\left(  \mathbf{r}\right)  =\theta\left(
\left|  \mathbf{r}\right|  -R_{i}\right)
\end{equation}%
\begin{equation}
\omega_{i}^{\left(  V2\right)  }\left(  \mathbf{r}\right)  =\frac{\mathbf{r}%
}{\left|  \mathbf{r}\right|  }\delta\left(  \left|  \mathbf{r}\right|
-R_{i}\right)
\end{equation}%
\begin{equation}
4\pi R_{i}^{2}\omega_{i}^{\left(  0\right)  }\left(  \mathbf{r}\right)  =4\pi
R_{i}\omega_{i}^{\left(  1\right)  }\left(  \mathbf{r}\right)  =\omega
_{i}^{\left(  2\right)  }\left(  \mathbf{r}\right)
\end{equation}%
\begin{equation}
4\pi R_{i}\omega_{i}^{\left(  V1\right)  }\left(  \mathbf{r}\right)
=\omega_{i}^{\left(  V2\right)  }\left(  \mathbf{r}\right)
\end{equation}
and $R_{i}$ is the radius of species $i$. $\delta$ is the Dirac delta function
and $\theta$ is the unit step function, $\theta\left(  x>0\right)  =0$ and
$\theta\left(  x\leq0\right)  =1$. Then the hard-sphere component $\mu
_{i}^{\text{HS}}$ of $\mu_{i}^{\text{ex}}$ is given by \cite{Rosenfeld89},
\cite{Rosenfeld93}%
\begin{equation}
\mu_{i}^{\text{HS}}\left(  \mathbf{x}\right)  =kT\sum_{\alpha}\int
\frac{\partial\Phi_{\text{HS}}}{\partial n_{\alpha}}\left(  \mathbf{x}%
^{\prime}\right)  \omega_{i}^{\left(  \alpha\right)  }\left(  \mathbf{x}%
-\mathbf{x}^{\prime}\right)  d\mathbf{x}^{\prime}. \label{mu_hs_eqn}%
\end{equation}

The integrals in Eqs. $\left(  \ref{n_def}\right)  $ and $\left(
\ref{mu_hs_eqn}\right)  $ simplify in the slab geometry. By analytically
integrating over two components of the three-dimensional vector $\mathbf{x}$
and assuming that the functions $c_{i}\left(  \mathbf{x}\right)  $ and
$\partial\Phi_{\text{HS}}/\partial n_{\alpha}\left(  \mathbf{x}\right)  $ are
constant over these two components, these integrals become \cite{Rosenfeld93}%
\begin{equation}
n_{\alpha}\left(  x\right)  =\sum_{i}\int_{x-R_{i}}^{x+R_{i}}c_{i}\left(
x^{\prime}\right)  W_{i}^{\left(  \alpha\right)  }\left(  x^{\prime}-x\right)
dx^{\prime} \label{pnl_def}%
\end{equation}%
\begin{equation}
\mu_{i}^{\text{HS}}\left(  x\right)  =kT\sum_{\alpha}\int_{x-R_{i}}^{x+R_{i}%
}\frac{\partial\Phi_{\text{HS}}}{\partial n_{\alpha}}\left(  x^{\prime
}\right)  W_{i}^{\left(  \alpha\right)  }\left(  x-x^{\prime}\right)
dx^{\prime} \label{mu_hs_1d}%
\end{equation}
where%
\begin{equation}
W_{i}^{\left(  2\right)  }\left(  r\right)  =2\pi R_{i}%
\end{equation}%
\begin{equation}
W_{i}^{\left(  3\right)  }\left(  r\right)  =\pi\left(  R_{i}^{2}%
-r^{2}\right)
\end{equation}%
\begin{equation}
W_{i}^{\left(  V2\right)  }\left(  r\right)  =2\pi r\left(  0,0,1\right)
\end{equation}%
\begin{equation}
4\pi R_{i}^{2}W_{i}^{\left(  0\right)  }\left(  r\right)  =4\pi R_{i}%
W_{i}^{\left(  1\right)  }\left(  r\right)  =W_{i}^{\left(  2\right)  }\left(
r\right)
\end{equation}%
\begin{equation}
4\pi R_{i}W_{i}^{\left(  V1\right)  }\left(  r\right)  =W_{i}^{\left(
V2\right)  }\left(  r\right)  .
\end{equation}

\subsection{Electrostatic component}

For the electrostatic component $\mu_{i}^{\text{ES}}$ of $\mu_{i}^{\text{ex}}$
we use a generalization of\ Rosenfeld's method \cite{Rosenfeld93} of expanding
$\Omega_{\text{ES}}$ in a functional Taylor series in powers of%
\begin{equation}
\Delta c_{i}\left(  \mathbf{x}\right)  =c_{i}\left(  \mathbf{x}\right)
-c_{i}^{\text{ref}}\left(  \mathbf{x}\right)
\end{equation}
where $c_{i}^{\text{ref}}\left(  \mathbf{x}\right)  $ is a given, nonconstant
reference concentration profile. (For simplicity, we choose this reference
fluid to be charge-neutral everywhere.) Then, up to second order,%
\begin{align}
\Omega_{\text{ES}}\left(  \left\{  c_{k}\left(  \mathbf{x}\right)  \right\}
\right)  \approx &  \Omega_{\text{ES}}\left(  \left\{  c_{k}^{\text{ref}%
}\left(  \mathbf{x}\right)  \right\}  \right)  +kT\sum_{i}\int\mu
_{i}^{\text{ES}}\left(  \left\{  c_{k}^{\text{ref}}\left(  \mathbf{x}\right)
\right\}  \right)  \Delta c_{i}\left(  \mathbf{x}\right)  d\mathbf{x}\\
&  -\frac{kT}{2}\sum_{i,j}\int\int c_{ij}^{\left(  2\right)  ,\text{ES}%
}\left(  \mathbf{x,x}^{\prime}\right)  \Delta c_{i}\left(  \mathbf{x}\right)
\Delta c_{j}\left(  \mathbf{x}^{\prime}\right)  d\mathbf{x}d\mathbf{x}%
^{\prime}\nonumber
\end{align}
where $c_{ij}^{\left(  2\right)  ,\text{ES}}$ is the electrostatic component
of the second-order direct correlation function of the reference fluid
\cite{Rosenfeld93}. Eq. $\left(  \ref{mu_es_from_gp}\right)  $ gives
\cite{Rosenfeld93}, \cite{Rosenfeld97}%
\begin{equation}
\mu_{i}^{\text{ES}}\left(  \mathbf{x}\right)  =\mu_{i}^{\text{ES}}\left(
\left\{  c_{k}^{\text{ref}}\left(  \mathbf{x}\right)  \right\}  \right)
-z_{i}e\phi\left(  \mathbf{x}\right)  -kT\sum_{j}\int c_{ij}^{\left(
2\right)  ,\text{ES}}\left(  \mathbf{x,x}^{\prime}\right)  \Delta c_{j}\left(
\mathbf{x}^{\prime}\right)  d\mathbf{x}^{\prime}. \label{mu_es_def}%
\end{equation}

\subsubsection{Non-reference fluid component}

To approximate the last term of Eq. $\left(  \ref{mu_es_def}\right)  $ we use
the method of Biben, Hansen, and Rosenfeld \cite{Rosenfeld97} which applies
the Mean Spherical Approximation (MSA) in the strong coupling limit.
Specifically, when two ions are far apart, they interact only through the
electric mean field and other correlations become insignificant. It is the MSA
to approximate all correlations outside of the contact radius%
\begin{equation}
R_{ij}=R_{i}+R_{j}%
\end{equation}
to be this asymptotic limit [Waisman\&Lebowitz]. In our case this gives%
\begin{equation}
kT\sum_{j}\int_{\left|  \mathbf{x-x}^{\prime}\right|  >R_{ij}}c_{ij}^{\left(
2\right)  ,\text{ES}}\left(  \mathbf{x,x}^{\prime}\right)  c_{j}\left(
\mathbf{x}^{\prime}\right)  d\mathbf{x}^{\prime}\approx-z_{i}e\phi\left(
\mathbf{x}\right)
\end{equation}
or \cite{Rosenfeld97}%
\begin{equation}
\mu_{i}^{\text{ES}}\left(  \mathbf{x}\right)  =\mu_{i}^{\text{ES}}\left(
\left\{  c_{k}^{\text{ref}}\left(  \mathbf{x}\right)  \right\}  \right)
-kT\sum_{j}\int_{\left|  \mathbf{x-x}^{\prime}\right|  \leq R_{ij}}%
c_{ij}^{\left(  2\right)  ,\text{ES}}\left(  \mathbf{x,x}^{\prime}\right)
\Delta c_{j}\left(  \mathbf{x}^{\prime}\right)  d\mathbf{x}^{\prime}.
\label{mu_es_3d}%
\end{equation}

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\pi R_{i}^{2}$ and $z_{j}e/4\pi R_{j}^{2}$, respectively
\cite{Rosenfeld93}, \cite{Rosenfeld97}:%
\begin{align}
&  kTc_{ij}^{\left(  2\right)  ,\text{ES}}\left(  \mathbf{x,x}^{\prime}\right)
\nonumber\\
&  =\frac{z_{i}z_{j}e^{2}}{4\pi\epsilon\left(  \mathbf{x}\right)  \epsilon
_{0}}\frac{1}{\left|  \mathbf{x-x}^{\prime}\right|  }-\frac{z_{i}z_{j}e^{2}%
}{8\pi\epsilon\left(  \mathbf{x}\right)  \epsilon_{0}}\frac{1}{\left|
\mathbf{x-x}^{\prime}\right|  }\left(  \left|  \mathbf{x-x}^{\prime}\right|
\frac{R_{i}+R_{j}}{R_{i}R_{j}}-\frac{\left|  \mathbf{x-x}^{\prime}\right|
^{2}}{2R_{i}R_{j}}-\frac{\left(  R_{i}-R_{j}\right)  ^{2}}{2R_{i}R_{j}}\right)
\\
&  =\frac{z_{i}z_{j}e^{2}}{16\pi\epsilon\left(  \mathbf{x}\right)
\epsilon_{0}}\frac{1}{R_{i}R_{j}}\frac{1}{\left|  \mathbf{x-x}^{\prime
}\right|  }\left[  \left|  \mathbf{x-x}^{\prime}\right|  -\left(  R_{i}%
+R_{j}\right)  \right]  ^{2} \label{cbar_eqn}%
\end{align}
where we have used the simplifying assumption that the dielectric coefficient
is approximately constant with a value of $\epsilon\left(  \mathbf{x}\right)
$. This is exact in the case of a constant dielectric coefficient $\epsilon$.
In the slab geometry, Eq. $\left(  \ref{mu_es_3d}\right)  $ becomes%
\begin{align}
\mu_{i}^{\text{ES}}\left(  x\right)  =  &  \mu_{i}^{\text{ES}}\left(  \left\{
c_{k}^{\text{ref}}\left(  x\right)  \right\}  \right)  -\sum_{j}%
\frac{z_{i}z_{j}e^{2}}{8\epsilon\left(  x\right)  \epsilon_{0}}\frac{1}%
{R_{i}R_{j}}\label{mu_es_1d}\\
&  \times\int_{x-R_{ij}}^{x+R_{ij}}c_{i}\left(  x^{\prime}\right)  \left(
\frac{1}{3}\sigma_{ij}^{3}+\sigma_{ij}\left(  x^{\prime}-x\right)
^{2}-\frac{1}{3}\left|  x^{\prime}-x\right|  ^{3}-\sigma_{ij}^{2}\left|
x^{\prime}-x\right|  \right)  dx^{\prime}.\nonumber
\end{align}

\subsubsection{Reference fluid component}

The only remaining component of the excess chemical potential is the
electrostatic excess chemical potential of the reference fluid, $\mu
_{i}^{\text{ES}}\left(  \left\{  c_{k}^{\text{ref}}\left(  x\right)  \right\}
\right)  $ in Eq. $\left(  \ref{mu_es_1d}\right)  $. In Rosenfeld's
implementation \cite{Rosenfeld93} 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.

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 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 Gauss-filtered
to produce continuous concentration profiles. A charge-neutral reference fluid
is produced by adjusting the anion concentrations.

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 \cite{NGHE}. The connecting regions are not homogeneous fluids and to
describe them correctly would require a 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.

\section{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 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 \cite{NGHE},
[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:%
\begin{equation}
\frac{d\mu_{j}}{dx}=0. \label{confined_equilib}%
\end{equation}
It is not a priori known, however, what the value of the chemical potential
$\mu_{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.
$\left(  \ref{confined_equilib}\right)  $, for each confined species we solve
the equation%
\begin{equation}
\int_{x_{L}}^{x_{R}}c_{j}\left(  x\right)  dx=\frac{N_{j}}{V}
\label{confined_number}%
\end{equation}
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.

Because the structural particles are ions, they contribute to the electric
field and therefore are also included in the Poisson equation $\left(
\ref{Poisson_eqn}\right)  $ that calculates the electrostatic potential
experienced by the permeating ions. However, because there is finite number of
confined particles, if the electrostatic potential from Eq. $\left(
\ref{Poisson_eqn}\right)  $ 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. 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 approach we choose to alleviate this
interaction is to 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 $\mu_{j}$, the electrostatic potential $\phi_{j}$ for confined
species $j$ is calculated from the modified Poisson equation%
\begin{equation}
-\frac{\epsilon}{A\left(  x\right)  }\frac{d}{dx}\left(  A\left(  x\right)
\frac{d\phi_{j}}{dx}\left(  x\right)  \right)  =e\sum_{i\neq j}z_{i}%
c_{i}\left(  x\right)  +ez_{j}\frac{N_{j}-1}{N_{j}}c_{j}\left(  x\right)
\label{confined_Poisson}%
\end{equation}
where the concentration of confined species $j$ is reduced by $1/N_{j}$ to
approximately eliminate the self-interaction.

\section{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. $\left(  \ref{NP_eqn}\right)  $, $\left(
\ref{continuity_eqn}\right)  $, $\left(  \ref{Poisson_eqn}\right)  $, $\left(
\ref{mu_hs_1d}\right)  $, $\left(  \ref{mu_es_1d}\right)  $, $\left(
\ref{confined_equilib}\right)  $, $\left(  \ref{confined_number}\right)  $,
and $\left(  \ref{confined_Poisson}\right)  $.

\subsection{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
they are baths 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\left(  x\right)  $, the area available for flux
\cite{capnp}.

To discretize the differential equations $\left(  \ref{NP_eqn}\right)  $,
$\left(  \ref{continuity_eqn}\right)  $, $\left(  \ref{Poisson_eqn}\right)  $,
and $\left(  \ref{confined_Poisson}\right)  $, these equations were rewritten
as follows:%
\begin{equation}
0=\frac{d}{dx}\left[  F_{i}\left(  x\right)  \frac{d}{dx}\left(
u_{i}+\frac{z_{i}e}{kT}\phi\left(  x\right)  +\frac{1}{kT}\mu_{i}^{0}\left(
x\right)  +\frac{1}{kT}\mu_{i}^{\text{ex}}\left(  x\right)  \right)  \right]
\label{predisc_NP}%
\end{equation}%
\begin{equation}
u_{i}\left(  x\right)  =\ln\left[  c_{i}\left(  x\right)  \right]
\end{equation}%
\begin{equation}
F_{i}\left(  x\right)  =\frac{D_{i}\left(  x\right)  A\left(  x\right)  }%
{\min_{x}\left\{  D_{i}\left(  x\right)  A\left(  x\right)  \right\}  }%
\exp\left[  u_{i}\left(  x\right)  \right]
\end{equation}%
\begin{equation}
-\frac{\epsilon\epsilon_{0}kT}{N_{A}e^{2}d^{2}}\frac{1}{A\left(  x\right)
}\frac{d}{dx}\left(  A\left(  x\right)  \frac{d\phi}{dx}\right)  =\sum
_{i}z_{i}\exp\left[  u_{i}\left(  x\right)  \right]  \label{predisc_Poisson}%
\end{equation}
where all concentrations have been scaled by Avogadro's number $N_{A}$ and all
lengths have been scaled by the system length $d$. (Eq. $\left(
\ref{confined_Poisson}\right)  $ may be rewritten in a similar form as Eq.
$\left(  \ref{predisc_Poisson}\right)  $).

Both the recast Nernst-Planck equation $\left(  \ref{predisc_NP}\right)  $ and
the recast Poisson equation $\left(  \ref{predisc_Poisson}\right)  $ have the
same form for the derivatives:%
\begin{equation}
\frac{d}{dx}\left(  f\left(  x\right)  \frac{dy}{dx}\right)  .
\end{equation}
We discretize this derivative on the interior points of the grid $\left\{
x_{0}=0,x_{1},x_{2},\ldots,x_{N-2},x_{N-1}=1\right\}  $ by
\begin{align}
&  \frac{d}{dx}\left(  f\left(  x\right)  \frac{dy}{dx}\right) \nonumber\\
&  \approx\frac{2}{h_{m}+h_{m+1}}\left[  f_{m+\frac{1}{2}}\frac{dy}{dx}\left(
x_{m+\frac{1}{2}}\right)  -f_{m-\frac{1}{2}}\frac{dy}{dx}\left(
x_{m-\frac{1}{2}}\right)  \right] \\
&  \approx\frac{2}{h_{m}+h_{m+1}}\left[  \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}}\right] \\
&  =\left(  f_{m}+f_{m+1}\right)  \left(  y_{m+1}-y_{m}\right)  \beta
_{1,m}-\left(  f_{m-1}+f_{m}\right)  \left(  y_{m}-y_{m-1}\right)
\beta_{-1,m}%
\end{align}
where%
\begin{equation}
f_{m}=f\left(  x_{m}\right)
\end{equation}%
\begin{equation}
y_{m}=y\left(  x_{m}\right)
\end{equation}%
\begin{equation}
\beta_{m,-1}=\frac{1}{h_{m}\left(  h_{m+1}+h_{m}\right)  }\qquad\beta
_{m,1}=\frac{1}{h_{m+1}\left(  h_{m+1}+h_{m}\right)  }%
\end{equation}%
\begin{equation}
h_{m}=x_{m}-x_{m-1}.
\end{equation}
This is finally rewritten as%
\begin{align}
&  \frac{d}{dx}\left(  f\left(  x\right)  \frac{dy}{dx}\right) \nonumber\\
&  \approx\beta_{-1,m}\left(  f_{m-1}+f_{m}\right)  y_{m-1}-\left[
\beta_{-1,m}\left(  f_{m-1}+f_{m}\right)  +\beta_{1,m}\left(  f_{m}%
+f_{m+1}\right)  \right]  y_{m}\\
&  +\beta_{1,m}\left(  f_{m}+f_{m+1}\right)  y_{m+1}.\nonumber
\end{align}
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},\ldots,x_{N-1}$.

The integrals in Eqs. $\left(  \ref{pnl_def}\right)  $, $\left(
\ref{mu_hs_1d}\right)  $, and $\left(  \ref{mu_es_1d}\right)  $ also 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 to 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. $\left(  \ref{pnl_def}\right)  $, $\left(
\ref{mu_hs_1d}\right)  $, and $\left(  \ref{mu_es_1d}\right)  $ one of the
functions can be written as a polynomial in the integration variable
$x^{\prime}$, each new integral in the sum can be evaluated analytically by
assuming that the other function in the integral ($c_{i}\left(  x^{\prime
}\right)  $ or $\partial\Phi_{\text{HS}}/\partial n_{\alpha}\left(  x^{\prime
}\right)  $) is a linear function in $x^{\prime}$ between consecutive grid nodes.

The remaining equations to be discretized are Eqs. $\left(
\ref{confined_equilib}\right)  $ and $\left(  \ref{confined_number}\right)  $
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. $\left(  \ref{confined_equilib}\right)  $ is
discretized as%
\begin{equation}
0=\mu_{j}\left(  x_{m}\right)  -\mu_{j}\left(  x_{m+1}\right)
\end{equation}
for $M_{L}\leq m\leq M_{R}-1$. For the last grid node of confinement
$x_{M_{R}}$, we apply Eq. $\left(  \ref{confined_number}\right)  $. 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.

\subsection{Solving the discretized system}

After discretizing the system of equations, the values of $u_{i}\left(
x_{m}\right)  =\ln\left[  c_{i}\left(  x_{m}\right)  \right]  $ for all
permeant species $i$ and $\phi\left(  x_{m}\right)  $ remain to be determined
at all interior grid points ($m=1,\ldots,N-2$), as well as $u_{j}\left(
x_{m}\right)  =\ln\left[  c_{i}\left(  x_{m}\right)  \right]  $ and $\phi
_{j}\left(  x_{m}\right)  $ for all confined species $j$ at all grid points to
which they are confined ($m=M_{L},\ldots,M_{R}$). We solved the system of
discretized equations using Newton's method \cite{Press} 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 \textit{all}
variables at \textit{every} iteration did not converge.

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. $\left(  \ref{confined_number}\right)  $ the band
becomes significantly more narrow. Because Eq. $\left(  \ref{confined_number}%
\right)  $ is a solution condition for Eq. $\left(  \ref{confined_equilib}%
\right)  $ 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. $\left(
\ref{confined_number}\right)  $ using an LU decomposition for band matrices
\cite{Press} and finally solving the complete matrix equation with the
Sherman-Morrison-Woodbury method \cite{Press}. 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}$.

\section{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.

\begin{thebibliography}{9}                                                                                                %

\bibitem {Hille}B.\ Hille, \textit{Ion Channels of Excitable Membranes,} 3rd
ed. (Sinauer Associates, Inc., Sunderland,\ MA, 2001).

\bibitem {Caccamo}C. Caccamo, J.-P. Hansen, and G. Stell, \textit{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).

\bibitem {Rosenfeld89}Y. Rosenfeld, Phys. Rev. Lett. \textbf{63}, 980 (1989).

\bibitem {Rosenfeld93}Y. Rosenfeld, J. Chem. Phys. \textbf{98}, 8126 (1993).

\bibitem {Rosenfeld97}Y. Rosenfeld, M.\ Schmidt, H. L\"{o}wen, and P.
Tarazona, Phys. Rev. E \textbf{55}, 4245 (1997).

\bibitem {Rosenfeld98}T. Biben, J.-P. Hansen, and Y. Rosenfeld, Phys. Rev. E
\textbf{57}, R3727 (1998).

\bibitem {SakmannNeher}B. Sakmann and E. Neher, eds., \textit{Single-Channel
Recording} (Plenum Press, New York, 1983).

\bibitem {originalSN}O. P. Hamill, A. Marty, E. Neher, B. Sakmann, and F. J.
Sigworth, Pfl\"{u}gers Arch. \textbf{391}, 85 (1981).

\bibitem {Romano}J. D. Romano and R. H. Price, Am. J. Phys. \textbf{64}, 1150 (1996).

\bibitem {Hansen}D. Goulding, S. Melchionna, and J.-P. Hansen, Phys. Chem.
Chem. Phys. \textbf{3}, 1644 (2001).

\bibitem {NGHE}W. Nonner, D. Gillespie, D. Henderson, and B. Eisenberg, J.
Phys. Chem. B \textbf{105}, 6427 (2001).

\bibitem {Chung}B. Corry, S. Keyucak, and S.-H. Chung, Biophys. J.
\textbf{78}, 2364 (2000).

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

\bibitem {frink}L. J. D. Frink, A. Thompson, and A. G. Salinger, J. Chem.
Phys. \textbf{112}, 7564 (2000).

\bibitem {thesis}D. Gillespie, Ph.D. thesis, Rush University, Chicago (1999).

\bibitem {capnp}W. Nonner and B. Eisenberg, Biophys. J. \textbf{75}, 1287 (1998).
\end{thebibliography}
\end{document}