%% This document created by Scientific Word (R) Version 3.5
\documentclass[12pt]{article}%
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{amsfonts}
\usepackage{amssymb}%
\setcounter{MaxMatrixCols}{30}
%TCIDATA{OutputFilter=latex2.dll}
%TCIDATA{Version=4.00.0.2309}
%TCIDATA{CSTFile=LaTeX article (bright).cst}
%TCIDATA{Created=Thursday, March 29, 2001 22:48:04}
%TCIDATA{LastRevised=Sunday, May 26, 2002 13:24:15}
%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}

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

\textbf{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. \textbf{XXXXX }In this paper we describe an engineering
approach to selective ion conduction through ion channels. Specifically, we
develop a one-dimensional model of\textbf{\ 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 \textbf{charged spheres} confined to the
channel. In this \textbf{one dimensional }model\textbf{ (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 \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}

\label{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}. \textbf{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 \textbf{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 \textbf{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 \textbf{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\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 \textbf{ }interface
\textbf{that maintained a fixed volume for the pore}, 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:
\textbf{In our model, the interface separates regions of different pressure.
T}he 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\textbf{ 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 \cite{NGHE}. Here we consider the simple
case that the protein maintains a constant average pore volume.
\textbf{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
\textbf{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.}

\textbf{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): \textbf{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.}

\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. \textbf{Comment:} \textbf{please explain. I do not understand and
perhaps the reader will have the same problem.}

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

\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 \textbf{under these conditions} slab
geometry\ \ \textbf{can be used }in mean-field theories like the one we
consider here \cite{Chung}.
\end{enumerate}

\section{Modeling ion flux}

\label{Modeling Ion Flux}We describe steady-state flux by the constituitive
law \textbf{(Ohm's or Fick's law for electrolytes)}%
\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}
\textbf{and the continuity 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, \textbf{PNP}) 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}

\label{DFT}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. \textbf{The deBroglie
wavelength is [physical phrase or clause] and is given by [insert equation].}

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 \textbf{Eisenberg, Klosek and Schuss,
REF who studied systems that might not be in local equilibrium,} and by Frink,
Thompson, and Salinger \cite{frink} who studied the steady-state diffusion of
colored hard spheres.

\subsection{Hard-sphere component}

\label{Hard-sphere component}For the hard-sphere component, we use the
\textquotedblleft antisymmetrized\textquotedblright\ excess free energy
density \cite[Eq. (27)]{Rosenfeld97} \textbf{Comment: The reader will want to
know why you chose this.}
\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\vert \mathbf{r}\right\vert -R_{i}\right)
\end{equation}%
\begin{equation}
\omega_{i}^{\left(  3\right)  }\left(  \mathbf{r}\right)  =\theta\left(
\left\vert \mathbf{r}\right\vert -R_{i}\right)
\end{equation}%
\begin{equation}
\omega_{i}^{\left(  V2\right)  }\left(  \mathbf{r}\right)  =\frac{\mathbf{r}%
}{\left\vert \mathbf{r}\right\vert }\delta\left(  \left\vert \mathbf{r}%
\right\vert -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 \ }

\label{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 \ }

\label{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
\textbf{in this limit}, they interact only through the electric mean field and
other correlations become insignificant.\textbf{The} MSA approximates all
correlations outside of the contact radius%
\begin{equation}
R_{ij}=R_{i}+R_{j}%
\end{equation}
\textbf{by } this asymptotic limit [Waisman\&Lebowitz]. In our case this gives%
\begin{equation}
kT\sum_{j}\int_{\left\vert \mathbf{x-x}^{\prime}\right\vert >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\vert \mathbf{x-x}^{\prime}\right\vert \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\vert \mathbf{x-x}^{\prime}\right\vert }-\frac{z_{i}%
z_{j}e^{2}}{8\pi\epsilon\left(  \mathbf{x}\right)  \epsilon_{0}}\frac
{1}{\left\vert \mathbf{x-x}^{\prime}\right\vert }\left(  \left\vert
\mathbf{x-x}^{\prime}\right\vert \frac{R_{i}+R_{j}}{R_{i}R_{j}}-\frac
{\left\vert \mathbf{x-x}^{\prime}\right\vert ^{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\vert \mathbf{x-x}^{\prime
}\right\vert }\left[  \left\vert \mathbf{x-x}^{\prime}\right\vert -\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$.
\textbf{Comment: the previous sentence needs to be explained. Why is it exact
then? }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\vert x^{\prime}-x\right\vert ^{3}-\sigma_{ij}^{2}\left\vert
x^{\prime}-x\right\vert \right)  dx^{\prime}.\nonumber
\end{align}


\subsubsection{Reference fluid component}

\label{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 \textquotedblleft intrinsic\textquotedblright\ reference fluid
compositions. This is the reason for choosing a non-constant reference fluid.

\textbf{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 \textbf{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 Gauss-filtered \textbf{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.

\textbf{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 \cite{NGHE}. The connecting regions are not homogeneous fluids and to
describe them correctly would require a (\textbf{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\textbf{, including its component the excess
electrostatic excess chemical potential.}

\section{Modeling the selectivity filter}

\label{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
\textquotedblleft EEEE\textquotedblright\ 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
\textquotedblleft DEKA\textquotedblright\ locus of the voltage-gated sodium
channel), then the channel becomes sodium selective [REFs].

In equilibrium studies, these \textquotedblleft structural
charges\textquotedblright\ have been modeled as \textbf{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
\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.
\textbf{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 $\left(
\ref{Poisson_eqn}\right)  $ that calculates the electrostatic potential
experienced by the permeating ions. However, because there is \textbf{a
}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. \textbf{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. \textbf{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 \textbf{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 $\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.

\textbf{Comment: this is a brilliant solution to the problem. Now let's see
Zeev prove it is right!!!}

\section{Numerical implementation}

\label{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}

\label{Discretization}The grid we chose was nonuniform because the baths are
necessarily\textbf{  }much longer than the channel to which the structural
charges are confined and, because \textbf{XXX} baths \textbf{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\left(
x\right)  $, the area available for flux \cite{capnp}.

\textbf{The differential equations XXX} $\left(  \ref{NP_eqn}\right)  $,
$\left(  \ref{continuity_eqn}\right)  $, $\left(  \ref{Poisson_eqn}\right)  $,
and $\left(  \ref{confined_Poisson}\right)  $, are rewritten as follows
\textbf{before they are discretized}:%
\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)  $ \textbf{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 \textbf{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. $\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}

\label{Solving the discretized system} After \textbf{XXX} the system of
equations \textbf{is discretized}, 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 \textbf{only converged (over the necessary domain of parameters) if
they update} \textit{all} variables at \textit{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. $\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}$.

\textbf{Comment: Please give an estimate of the time the calculation takes.}

\section{Example: a calcium channel}

\label{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}{99}                                                                                               %


\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}