%% 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=Tuesday, May 28, 2002 13:26:07}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{<META NAME="DocumentShell" CONTENT="General\SW\Blank - Standard LaTeX Article">}
%TCIDATA{ComputeDefs=
%$f\left(  r\right)  =-\frac{z_{i}e}{r^{3}}$
%}
\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 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 \cite{Hille}.

Ion channels are proteins with functional groups that often are the side
chains of charged amino acids. Enough data on these functional groups and even
crystallographic channel structures are now available to apply specific
theories of permeation and selectivity. With such information and theories it
is possible to reverse engineer the devices that biologists observe, not only
with the goal of understanding their natural function, but also of controlling
their function beyond biology. In this paper we describe an engineering
approach to selective ion conduction through ion channels. Specifically, we
develop a one-dimensional model of the movement of ions between two baths of
fixed ionic concentrations and applied electrostatic potentials. The baths are
connected by a single\ open ion channel whose functional groups are
represented as ions confined to the channel. In this model (unlike some we
have studied \cite{LS}, \cite{CRC2}, \cite{capnp}) 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}

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

The geometry of two baths connected by a single ion channel embedded in a
membrane is shown in Fig. 1. This is the geometry of the most direct
experiments \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 controlled in
the channel, the resolution must be high there. However, the boundary
conditions of the problem must be applied far away from the channel to model
the problem correctly. This problem is three dimensional, but a reasonable
one-dimensional approximation can be developed that still includes the effects
of the finite volume of ions and functional groups of the channel protein. In
this section we derive an approximation that is fine grained in the
flux-limiting region (the channel) and coarse grained in the baths. This
approximation allows the problem to be solved on a desktop computer.

Starting in the baths, consider spherical shells that terminate at the lipid
and are perpendicular to the lipid as well as the $x$-axis (Fig. 1). The flux
of particles is perpendicular to these surfaces and therefore flux densities
scale inversely with the area $A\left(  x\right)  $ of these shells. (The
shells are indexed by their intersection with the $x$-axis.) This spherical
geometry is exact for each bath alone: ions traveling to and from the channel
must converge to a hemisphere given by the capture radius of the channel. As
we move from the baths into the channel, we continue to assume that particles
travel perpendicular to these spherical shells. Below we will argue that this
is an appropriate approximation inside the channel where the spherical shells
are flat, cross-sectional disks. In the transition regions this approximation
is also valid. For example, in the calculation of the electric field in these
transition regions, use of the spherical shells is equivalent to approximating
a flat-bottomed cone by a spherical cone \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 smooth, hard wall at the pore/protein interface that
maintained a fixed volume for the pore, then large radial packing effects
would result in the pore \cite{Hansen}. However, real proteins are not smooth
or rigid, and therefore we take the opposite view of a hard-wall interface:
the atoms of the channel protein are another liquid, forming a statistical
interface with the contents of the pore and providing a confining normal
pressure on the pore contents. Such a boundary is not described by a geometric
surface at a specific location, but rather by an \textit{average}
cross-sectional area or pore volume that may be constant or \textbf{vary with
conditions }depending on the compressibility of the protein \cite{NGHE}. Here
we consider the simple case that the protein maintains a constant
\textit{average} cross-section. \textbf{Comment: the word `variable' might
imply `vary with time' and that would be confusing. I find the word `variable'
often causes this problem for readers.}

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

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

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

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

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

\textbf{"Hard sphere interactions extend through the pore because....."}

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

\section{Modeling ion flux}

We describe steady-state flux by the constitutive law%
\begin{equation}
-J_{i}=\frac{1}{kT}D_{i}\left(  x\right)  A\left(  x\right)  \rho_{i}\left(
x\right)  \frac{d\mu_{i}}{dx} \label{NP_eqn_alt}%
\end{equation}
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}\left(  x\right)  $ 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. This
model was recently derived from a Poisson-Langevin system \cite{SNE}.

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)  $ from particle interactions, 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[  \frac{\rho_{i}\left(  x\right)  }{\rho_{\text{scale}}}\right]
\label{ideal_def}%
\end{equation}
where $\phi\left(  x\right)  $ is the local electrostatic potential, $z_{i}$
is the valence of species $i$, and $\rho_{\text{scale}}$ is some
characteristic number density. $\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, PNP) system:
\begin{equation}
-J_{i}=D_{i}\left(  x\right)  A\left(  x\right)  \left[  \frac{d\rho_{i}}%
{dx}+\frac{1}{kT}\rho_{i}\left(  x\right)  \left(  z_{i}e\frac{d\phi}%
{dx}+\frac{d\mu_{i}^{0}}{dx}+\frac{d\mu_{i}^{\text{ex}}}{dx}\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}\right)  =e\sum_{i}z_{i}\rho_{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}\right)  =e\sum_{i}z_{i}\rho_{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\{  \rho_{k}\left(  \mathbf{x}\right)  \right\}  \right)
=\Omega_{\text{id}}\left(  \left\{  \rho_{k}\left(  \mathbf{x}\right)
\right\}  \right)  +\Omega_{\text{HS}}\left(  \left\{  \rho_{k}\left(
\mathbf{x}\right)  \right\}  \right)  +\Omega_{\text{ES}}\left(  \left\{
\rho_{k}\left(  \mathbf{x}\right)  \right\}  \right)  \label{grand_potl}%
\end{equation}
with%
\begin{equation}
\Omega_{\text{id}}\left(  \left\{  \rho_{k}\left(  \mathbf{x}\right)
\right\}  \right)  =kT\sum_{i}\int\rho_{i}\left(  \mathbf{x}\right)  \left[
\ln\left(  \frac{\rho_{i}\left(  \mathbf{x}\right)  }{\rho_{\text{scale}}%
}\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}
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\{  \rho_{k}\left(  \mathbf{x}^{\prime}\right)  \right\}
\right)  }{\delta\rho_{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\{  \rho_{k}\left(  \mathbf{x}^{\prime}\right)  \right\}
\right)  }{\delta\rho_{i}\left(  \mathbf{x}\right)  }. \label{mu_es_from_gp}%
\end{equation}


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

\subsection{Hard-sphere component}

For the hard-sphere component, we use the ``antisymmetrized'' excess free
energy density \cite[Eq. (27)]{Rosenfeld97} because it is the best currently
available for particles in confined geometries:
\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\rho_{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)  \label{vol_weight}%
\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 $\rho_{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}}\rho_{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\rho_{i}\left(  \mathbf{x}\right)  =\rho_{i}\left(  \mathbf{x}\right)
-\rho_{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. Then, up to second order,%
\begin{align}
\Omega_{\text{ES}}\left(  \left\{  \rho_{k}\left(  \mathbf{x}\right)
\right\}  \right)  \approx &  \Omega_{\text{ES}}\left(  \left\{  \rho
_{k}^{\text{ref}}\left(  \mathbf{x}\right)  \right\}  \right)  +kT\sum_{i}%
\int\mu_{i}^{\text{ES}}\left(  \left\{  \rho_{k}^{\text{ref}}\left(
\mathbf{x}\right)  \right\}  \right)  \Delta\rho_{i}\left(  \mathbf{x}\right)
d\mathbf{x}\label{ES_expand}\\
&  -\frac{kT}{2}\sum_{i,j}\int\int c_{ij}^{\left(  2\right)  ,\text{ES}%
}\left(  \mathbf{x,x}^{\prime}\right)  \Delta\rho_{i}\left(  \mathbf{x}%
\right)  \Delta\rho_{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\{  \rho_{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\rho
_{j}\left(  \mathbf{x}^{\prime}\right)  d\mathbf{x}^{\prime}.
\label{mu_es_def}%
\end{equation}


\subsubsection{Reference fluid component}

\label{reference_sec}

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 regions of
the system are very diverse. Channels are embedded in a membrane that
separates two bulk fluids whose ionic concentrations are different,
experimentally controlled, and can range from $<1$ $\mu$M to 2 M. More
importantly, inside the channel there is a very high concentration of charges
from the amino acid residues of the channel protein. These attract equally
high concentrations of counter charge with the total concentration of tens of
molar. Thus, each section of the geometry has a different ``intrinsic''
reference fluid composition. This is the reason for choosing a non-uniform
reference fluid.

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

To determine the $\rho_{k}^{\text{ref}}\left(  \mathbf{x}\right)  $, we start
by defining%
\begin{equation}
\bar{\rho}_{k}\left(  \mathbf{x}\right)  =A\rho_{k}^{\ast}\left(
\mathbf{x}\right)  \qquad\left(  z_{k}\geq0\right)
\end{equation}
and%
\begin{equation}
\bar{\rho}_{k}\left(  \mathbf{x}\right)  =AB\rho_{k}^{\ast}\left(
\mathbf{x}\right)  \qquad\left(  z_{k}<0\right)
\end{equation}
with the constants $A$ and $B$ determined by requiring the $\left\{  \bar
{\rho}_{k}\left(  \mathbf{x}\right)  \right\}  $ to be charge-neutral and to
have the same ionic strength at each location $\mathbf{x}$ as the $\left\{
\rho_{k}^{\ast}\left(  \mathbf{x}\right)  \right\}  $; that is, we require
that%
\begin{equation}
0=\sum_{k}z_{k}\bar{\rho}_{k}\left(  \mathbf{x}\right)
\label{ref_neutral_cond}%
\end{equation}
and%
\begin{equation}
\sum_{k}z_{k}^{2}\bar{\rho}_{k}\left(  \mathbf{x}\right)  =\sum_{k}z_{k}%
^{2}\rho_{k}^{\ast}\left(  \mathbf{x}\right)  .\label{ref_is_cond}%
\end{equation}
These conditions give \textbf{XXXX}%
\begin{equation}
A=\frac{\sum_{k}z_{k}^{2}\rho_{k}^{\ast}\left(  \mathbf{x}\right)  }%
{\sum_{z_{k}\geq0}z_{k}^{2}\rho_{k}^{\ast}\left(  \mathbf{x}\right)
+B\sum_{z_{k}<0}z_{k}^{2}\rho_{k}^{\ast}\left(  \mathbf{x}\right)  }%
\end{equation}
and%
\begin{equation}
B=\frac{\sum_{z_{k}\geq0}z_{k}\rho_{k}^{\ast}\left(  \mathbf{x}\right)  }%
{\sum_{z_{k}<0}\left\vert z_{k}\right\vert \rho_{k}^{\ast}\left(
\mathbf{x}\right)  }.
\end{equation}


The $\rho_{k}^{\text{ref}}\left(  \mathbf{x}\right)  $ are then defined as
spatially averaged $\bar{\rho}_{k}\left(  \mathbf{x}\right)  $:%
\begin{equation}
\rho_{k}^{\text{ref}}\left(  \mathbf{x}\right)  =\int\bar{\rho}_{k}\left(
\mathbf{x}^{\prime}\right)  w\left(  \mathbf{x}^{\prime},\mathbf{x}\right)
d\mathbf{x}^{\prime}%
\end{equation}
where the weight function $w$ is a normalized form of the volume weight
$\omega_{i}^{\left(  3\right)  }$ (Eq. $\left(  \ref{vol_weight}\right)  $):
\begin{equation}
w\left(  \mathbf{x}^{\prime},\mathbf{x}\right)  =\frac{\theta\left(  \left|
\mathbf{x}^{\prime}-\mathbf{x}\right|  -R_{\text{filter}}\left(
\mathbf{x}\right)  \right)  }{\frac{4\pi}{3}R_{\text{filter}}^{3}\left(
\mathbf{x}\right)  }.
\end{equation}
The radius of the sphere over which we average is an approximation of the
electrostatic length scale, namely the capacitance radius (that is, the ion
radius plus the screening length):%
\begin{equation}
R_{\text{filter}}\left(  \mathbf{x}\right)  =s^{\ast}\left(  \mathbf{x}%
\right)  +\frac{\sum_{k}\bar{\rho}_{k}\left(  \mathbf{x}\right)  R_{k}}%
{\sum_{k}\bar{\rho}_{k}\left(  \mathbf{x}\right)  }.
\end{equation}
The averaging produces a kind of non-local reference concentration, removing
discontinuities that may have been in the $\left\{  \bar{\rho}_{k}\left(
\mathbf{x}\right)  \right\}  $. Futhermore, because we use the same filter and
filter length for each species $k$, the resulting $\left\{  \rho
_{k}^{\text{ref}}\left(  \mathbf{x}\right)  \right\}  $ are charge-neutral at
every point since the $\left\{  \bar{\rho}_{k}\left(  \mathbf{x}\right)
\right\}  $ are.

To calculate the ES excess chemical potential of the reference fluid, we apply
the MSA for bulk fluids at every point [REFs]. This not only results in values
for $\mu_{i}^{\text{ES}}\left(  \left\{  \rho_{k}^{\text{ref}}\left(
\mathbf{x}\right)  \right\}  \right)  $, but also for the screening length,
which we now redefine as%
\begin{equation}
s\left(  \mathbf{x}\right)  =\frac{1}{2\Gamma\left(  \mathbf{x}\right)  }
\label{new_screen_eq}%
\end{equation}
where $\Gamma\left(  \mathbf{x}\right)  $ is the MSA parameter $\Gamma$ that
has been calculated at every point $\mathbf{x}$. This then defines the
capacitance length of each ion species:%
\begin{equation}
\lambda_{k}\left(  \mathbf{x}\right)  =R_{k}+\lambda\left(  \mathbf{x}\right)
. \label{cap_len_eq}%
\end{equation}


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

\subsubsection{Non-reference fluid component}

Calculating the last term of Eq. $\left(  \ref{mu_es_def}\right)  $ requires
the correlation function of the reference fluid $c_{ij}^{\left(  2\right)
,\text{ES}}\left(  \mathbf{x,x}^{\prime}\right)  $. We start by approximating
the direct correlation function with the MSA \cite{Waisman&Lebowitz}:%
\begin{equation}
c_{ij}^{\left(  2\right)  ,\text{ES}}\left(  \mathbf{x,x}^{\prime}\right)
=-\frac{1}{kT}\psi_{ij}\left(  \mathbf{x,x}^{\prime}\right)
\end{equation}
if%
\begin{equation}
\left|  \mathbf{x-x}^{\prime}\right|  >R_{ij}=R_{i}+R_{j}%
\end{equation}
where $\psi_{ij}\left(  \mathbf{x,x}^{\prime}\right)  $ is the interaction
potential of two point particles with charge $z_{i}e$ and $z_{j}e$ located at
$\mathbf{x}$ and $\mathbf{x}^{\prime}$, respectively. In the case of uniform
dielectric coefficient,%
\begin{equation}
\psi_{ij}\left(  \mathbf{x,x}^{\prime}\right)  =\frac{z_{i}z_{j}e^{2}}%
{4\pi\epsilon\epsilon_{0}}\frac{1}{\left|  \mathbf{x-x}^{\prime}\right|  }.
\end{equation}
Since%
\begin{equation}
kT\sum_{j}\int\psi_{ij}\left(  \mathbf{x,x}^{\prime}\right)  \rho_{j}\left(
\mathbf{x}^{\prime}\right)  d\mathbf{x}^{\prime}=z_{i}e\phi\left(
\mathbf{x}\right)  ,
\end{equation}
we have%
\begin{align}
\mu_{i}^{\text{ES}}\left(  \mathbf{x}\right)   &  =\mu_{i}^{\text{ES}}\left(
\left\{  \rho_{k}^{\text{ref}}\left(  \mathbf{x}\right)  \right\}  \right)
\label{mu_es_3d}\\
&  -kT\sum_{j}\int_{\left|  \mathbf{x-x}^{\prime}\right|  \leq R_{ij}}\left(
c_{ij}^{\left(  2\right)  ,\text{ES}}\left(  \mathbf{x,x}^{\prime}\right)
+\psi_{ij}\left(  \mathbf{x,x}^{\prime}\right)  \right)  \Delta\rho_{j}\left(
\mathbf{x}^{\prime}\right)  d\mathbf{x}^{\prime}\nonumber
\end{align}
where the core-overlap electrostatic correlations remain to be determined. A
reasonable choice is the bulk MSA ES correlation function which can be
expressed analytically even for mixtures \cite{Hiroike}, \cite{BlumRosenfeld}.
We choose \textit{not} to use this formulation because the bulk MSA ES
correlation function involves a uniform dielectric coefficient. While we use a
uniform dielectric in this paper, to study ion permeation through ion channels
more correctly, one must eventually use an inhomogeneous dielectric coeffient
\cite{NGHE}. We will do this in future work and this will require an
appropriate ES correlation function. The approximate ES correlation function
we use in this paper is one that can be extended to the inhomogeneous
dielectric coefficient case; we test it here for the uniform dielectric
coefficient case, before applying it to more general situations. Specifically,
we use the approximation given by Blum and Rosenfeld \cite[Eq. (56)]%
{BlumRosenfeld} that the ES correlation function $c_{ij}^{\left(  2\right)
,\text{ES}}$ is the negative of the interaction potential of two spherical
shells $\Theta_{i}$ and $\Theta_{j}$ with each $\Theta_{k}$ having charge
equal and opposite of the ion species $k$ and radius $\lambda_{k}$, the
capacitance radius of ion species $k$ defined by Eq. $\left(  \ref{cap_len_eq}%
\right)  $.

With this ansatz, the short-range electrostatic correlations reduce to the
interaction energy of the corresponding point charges $z_{i}e$ and $z_{j}e$
minus the interaction energy of particles of radius $R_{i}$ and $R_{j}$ with
surface charge $z_{i}e/4\pi R_{i}^{2}$ and $z_{j}e/4\pi R_{j}^{2}$,
respectively \cite{Rosenfeld93}, \cite{Rosenfeld97}:%
\begin{equation}
kTc_{ij}^{\left(  2\right)  ,\text{ES}}\left(  \mathbf{x,x}^{\prime}\right)
=-\frac{z_{i}z_{j}e^{2}}{8\pi\epsilon\epsilon_{0}}\frac{1}{\left|
\mathbf{x-x}^{\prime}\right|  }\left(  \left|  \mathbf{x-x}^{\prime}\right|
\frac{\lambda_{i}+\lambda_{j}}{\lambda_{i}\lambda_{j}}-\frac{\left|
\mathbf{x-x}^{\prime}\right|  ^{2}}{2\lambda_{i}\lambda_{j}}-\frac{\left(
\lambda_{i}-\lambda_{j}\right)  ^{2}}{2\lambda_{i}\lambda_{j}}\right)
\label{dcf_sr_eqn}%
\end{equation}
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\{
\rho_{k}^{\text{ref}}\left(  x\right)  \right\}  \right)  -\sum_{j}\frac
{z_{i}z_{j}e^{2}}{8\epsilon\epsilon_{0}}\frac{1}{\lambda_{i}\lambda_{j}%
}\label{mu_es_1d}\\
&  \times\int_{x-R_{ij}}^{x+R_{ij}}c_{i}\left(  x^{\prime}\right)  \left(
\frac{1}{3}R_{ij}^{3}+\lambda_{ij}^{2}R_{ij}-\lambda_{ij}\allowbreak
R_{ij}^{2}+\lambda_{ij}(x^{\prime}-x)^{2}\right. \nonumber\\
&  \left.  -\frac{1}{3}\left|  x^{\prime}-x\right|  ^{3}-\lambda_{ij}%
^{2}\allowbreak\left|  x^{\prime}-x\right|  \right)  dx^{\prime}\nonumber
\end{align}
where%
\begin{equation}
\lambda_{ij}=\lambda_{i}+\lambda_{j}.
\end{equation}


\subsection{Iterating to self-consistency}

\label{iteration_sec}

\textbf{I think this section is fabulous science and very good writing.}

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

\subsection{Tests of the solution method}

\label{calibrations_sec}

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

For these \textbf{test }simulations, the equations were solved as described
below (Sec. \ref{numerics_sec}) except that because\ \textbf{the test XXX} is
an equilibrium problem we did not solve the Nernst-Planck equations $\left(
\ref{NP_eqn}\right)  $, but rather%
\begin{equation}
\mu_{k}\left(  x\right)  =\mu_{k}%
\end{equation}
for ion species where the total electrochemical potential $\mu_{k}$ was given
for each ion. The results of the simulations (shown in Figs. 2-4) are quite good.

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

\section{Modeling the selectivity filter}

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

In equilibrium studies, these essential \textquotedblleft structural
charges\textquotedblright\ have been modeled as ions that are allowed to move
freely inside the selectivity filter, but are \textbf{confined there. They 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 what the value of the chemical potential $\mu_{j}$
is. What is known, however, is that the number of particles of each confined
species $j$ is fixed ($N_{j}$). Therefore, in conjunction with Eq. $\left(
\ref{confined_equilib}\right)  $, for each confined species we solve the
equation%
\begin{equation}
\int_{x_{L}}^{x_{R}}\rho_{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{COMMENT: I continue to think something should be said about this
equation differing from the constant volume treatment of Hansen. Please
explain to me if I am wrong, why I am wrong. If I am right, i.e., that this is
very different in the two treatments, a remark should be made however
politely. COMMENT: for next paragraph. I liked what I suggested better. I do
not mind being over-ruled but please explain why. I give what I wrote again
for this section below, separated by marker lines}

Structural particles that are ions contribute to the electric field and
therefore are also included in the Poisson equation $\left(  \ref{Poisson_eqn}%
\right)  $ that calculates the electrostatic potential experienced by the
permeating ions. However, because there is a small number of confined
particles, if the electrostatic potential from Eq. $\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 includes a contribution from the particle itself. This is a problem
inherent to the mean-field approach, but one that is significant only when
there are a small number \textbf{OF} species. \textbf{COMMENT: I still do not
understand the previous sentence except for the obvious that with an infinite
number of particles NOT species, we approach the mean field. I\ think it
essential to explain or reviewers will complain.}

We alleviate this self-interaction by calculating the electrochemical
potential of the confined particles with a different electrostatic potential
for each confined species that applies only to that confined particle species.
Specifically, in the calculation of $\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_{k\neq j}z_{k}\rho
_{k}\left(  x\right)  +ez_{j}\frac{N_{j}-1}{N_{j}}\rho_{j}\left(  x\right)
\label{confined_Poisson}%
\end{equation}
where the concentration of confined species $j$ is reduced by the factor
$1/N_{j}$ to approximately eliminate the self-interaction.

\textbf{\#\#\#\#\#\#\#\# PREVIOUS SUGGESTION FOLLOWS \#\#\#\#\#\#}

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

\textbf{......}%

\begin{equation}
\int_{x_{L}}^{x_{R}}\rho_{j}\left(  x\right)  dx=\frac{N_{j}}{V}%
\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)
\end{equation}
where the concentration of confined species $j$ is reduced by $1/N_{j}$ to
approximately eliminate theself-interaction.

\textbf{\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\# END OF PREVIOUS SUGGESTION
\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#}

\section{Numerical implementation}

\label{numerics_sec}

The inputs to the problem are the left and right far-bath concentrations of
the permeating particles (ions and solvent), the applied electrostatic
potential, and the number of structural (confined) ions, as well as the radii
and diffusion coeffients of all particles and the dielectric coefficient. The
outputs are the concentrations of all particles (permeating and structural)
and the electrostatic potential, both as functions of location. These are
found by simultaneously solving Eqs. $\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)  $. After these have been solved,
fluxes are computed from an integrated Nernst-Planck equation $\left(
\ref{NP_eqn}\right)  $:%
\begin{equation}
J_{i}=\frac{\rho_{i}\left(  L\right)  \exp\left(  \frac{1}{kT}\left[  \mu
_{i}\left(  L\right)  -\ln\left(  \rho_{i}\left(  L\right)  \right)  \right]
\right)  -\rho_{i}\left(  R\right)  \exp\left(  \frac{1}{kT}\left[  \mu
_{i}\left(  R\right)  -\ln\left(  \rho_{i}\left(  R\right)  \right)  \right]
\right)  }{\int\frac{\exp\left(  \frac{1}{kT}\left[  \mu_{i}\left(  x\right)
-\ln\left(  \rho_{i}\left(  x\right)  \right)  \right]  \right)  }%
{D_{i}\left(  x\right)  A\left(  x\right)  }dx}%
\end{equation}
where the integral is over the entire system and $\rho_{i}\left(  L/R\right)
$ are the given left ($L$) and right ($R$) bath concentrations (similarly for
the chemical potentials $\mu_{i}\left(  L/R\right)  $).

\subsection{Discretization}

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

The differential equations $\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 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[  \rho_{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)  $ 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. $\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 ($\rho_{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 the one-half rule. \textbf{COMMENT:
readers may wonder, as I do, what is the "one-half rule"??}

\subsection{Solving the discretized system}

After the system of equations is discretized, the values of $u_{i}\left(
x_{m}\right)  =\ln\left[  \rho_{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[  \rho_{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 \textbf{converge, COMMENT: comma
was inserted} or required a very large number of iterations. Using Newton's
method with a poor initial guess, 30 iterations are typical with a total
calculation time of one to two minutes, depending on the number ion species.
Many fewer iterations (4-10) are necessary when a better inital guess is
available, for example when calculating a current/voltage curve. Then a good
intial guess is a previous solution of the system with a different applied
voltage (up to approximately $0.1$ V difference).

Because of the specific functional dependence of the equations, it is possible
to write the Jacobian matrix as a band matrix. Furthermore, without the
discretized version of Eq. $\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}

\label{example_sec}

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

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

\section*{Acknowledgments}

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

\begin{thebibliography}{99}                                                                                               %


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

\bibitem {LS}\noindent D. P. Chen, J. Lear, and B. Eisenberg, Biophys. J.
\textbf{72}, 97 (1997).

\bibitem {CRC2}D. P. Chen, \textit{et al.,} Biophys.\ J. \textbf{76}, 1346 (1999).

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

\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 {SNE}Z. Schuss, B. Nadler, and B.\ Eisenberg, Phys. Rev. E
\textbf{64}, 036116 (2001).

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

\bibitem {EKS}R. S. Eisenberg, M. M. Klosek, and Z. Schuss, J. Chem. Phys.
\textbf{102}, 1767 (1995).

\bibitem {Waisman&Lebowitz}E. Waisman and J. L. Lebowitz, J. Chem. Phys.
\textbf{56}, 3086 (1972).

\bibitem {Hiroike}K. Hiroike, Mol. Phys. \textbf{33}, 1195 (1977).

\bibitem {BlumRosenfeld}L. Blum and Y. Rosenfeld, J. Stat. Phys. \textbf{63},
1177 (1991).

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

\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).
\end{thebibliography}


\end{document}