Modeling Charge-Sign Asymmetric Solvation Free Energies With Nonlinear Boundary Conditions
Jaydeep P. Bardhan
Matthew G. Knepley
We show that charge-sign-dependent asymmetric hydration can be modeled accurately using linear Poisson theory but replacing the standard electric-displacement boundary condition with a simple nonlinear boundary condition. Using a single multiplicative scaling factor to determine atomic radii from molecular dynamics Lennard-Jones parameters, the new model accurately reproduces MD free-energy calculations of hydration asymmetries for (i) monatomic ions, (ii) titratable amino acids in both their protonated and unprotonated states, and (iii) the Mobley ``bracelet'' and ``rod'' test problems [J. Phys. Chem. B, v. 112:2408, 2008]. Remarkably, the model also justifies the use of linear response expressions for charging free energies. Our boundary-element method implementation demonstrates the ease with which other continuum-electrostatic solvers can be extended to include asymmetry.
draft.tex
\documentclass[reprint, letterpaper, nobibnotes, aps, superscriptaddress,prb]{revtex4-1}

\usepackage{etoolbox}
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{color}
\usepackage{url}
\usepackage{multirow}
\usepackage{graphicx}        % standard LaTeX graphics tool
                             % when including figure files
%\usepackage{minted}



% good reviewers
%  Mike Gilson
%  Chris Fennell
%  Marilyn Gunner
%  Nathan Baker
%  Walter Rocchia
%  Emil Alexov
%  Stefan Boresch
%  Paul Tavan
%  Greg Schenter, PNNL
%  Jay Ponder
%  Shekhar Garde
%  Monte Pettitt
%  

\def\vf{{\bf f}}
\def\vrp{{\bf r}'}
\def\vr{{\bf r}}
\def\vu{{\bf u}}
\def\vv{{\bf v}}
\def\vx{{\bf x}}
\def\vn{{\bf n}}
\def\code#1{{\tt #1}}
\def\enum#1{{\sc #1}}
\def\class#1{{\tt\bf #1}}
\def\function#1{\code{#1}}

% Color macros
\newcommand{\blue}{\textcolor{blue}}
\newcommand{\green}{\textcolor{green}}
\newcommand{\red}{\textcolor{red}}
\newcommand{\brown}{\textcolor{brown}}
\newcommand{\cyan}{\textcolor{cyan}}
\newcommand{\magenta}{\textcolor{magenta}}
\newcommand{\yellow}{\textcolor{yellow}}

% Author notes
%\usepackage[colorinlistoftodos, textwidth=4cm, shadow]{todonotes}
%\newcommand{\note}[1]{\todo[color=orange!40,inline]{{\bf Note:} #1}}
%\newcommand{\jay}[1]{\todo[color=blue!40,inline]{{\bf Jay:} #1}}
%\newcommand{\matt}[1]{\todo[color=green!40,inline]{{\bf Matt:} #1}}
%\newcommand{\fixme}[1]{\todo[color=red!40,inline]{{\bf FIXME:} #1}}
% \newcommand{\note}[1]{}
% \newcommand{\andy}[1]{}
% \newcommand{\matt}[1]{}
% \newcommand{\fixme}[1]{}
\let\comment = \note

\newcommand{\arxiv}{\let\arxivFig\matt}

\def\Xint#1{\mathchoice
   {\XXint\displaystyle\textstyle{#1}}%
   {\XXint\textstyle\scriptstyle{#1}}%
   {\XXint\scriptstyle\scriptscriptstyle{#1}}%
   {\XXint\scriptscriptstyle\scriptscriptstyle{#1}}%
   \!\int}
\def\XXint#1#2#3{{\setbox0=\hbox{$#1{#2#3}{\int}$}
     \vcenter{\hbox{$#2#3$}}\kern-.5\wd0}}
\def\ddashint{\Xint=}
\def\dashint{\Xint-}

\newtoggle{fulltitlepage}
\settoggle{fulltitlepage}{true}

\begin{document}

\title{Modeling Charge-Sign Asymmetric Solvation Free Energies With Nonlinear Boundary Conditions}

\iftoggle{fulltitlepage}{
\author{Jaydeep P. Bardhan}
\affiliation{Dept.~of Electrical and Computer Enginering, Northeastern University, Boston MA 02115}
\author{Matthew G. Knepley}
\affiliation{Computation Institute, The University of Chicago, Chicago IL 60637}
\begin{abstract}
We show that charge-sign-dependent asymmetric hydration can be modeled
accurately using linear Poisson theory but replacing the standard
electric-displacement boundary condition with a simple nonlinear
boundary condition.  Using a single multiplicative scaling factor to
determine atomic radii from molecular dynamics Lennard-Jones
parameters, the new model accurately reproduces MD free-energy
calculations of hydration asymmetries for (i) monatomic ions, (ii)
titratable amino acids in both their protonated and unprotonated
states, and (iii) the Mobley ``bracelet'' and ``rod'' test problems
[J. Phys. Chem. B, v. 112:2408, 2008].  Remarkably, the model also
justifies the use of linear response expressions for charging free
energies.  Our boundary-element method implementation demonstrates the
ease with which other continuum-electrostatic solvers can be extended
to include asymmetry.
\end{abstract}
}{
}

\maketitle

\section{Introduction}\label{sec:intro}
Implicit-solvent models represent an intuitive and fast approach to
understand molecular solvation~\cite{Roux99,Sharp90,Davis90,Tomasi94}, and
have a rigorous statistical-mechanical interpretation as an
approximation to the potential of mean force (PMF) experienced by a
molecular solute due to the surrounding solvent
molecules~\cite{Roux99}.  The PMF is usually decomposed into non-polar
and electrostatic terms, the latter of which are often modeled using
macroscopic continuum models based on the Poisson--Boltzmann
partial-differential equation (PDE).  Continuum models approximate the
free energy required to grow the solute charge distribution into the
solute
cavity~\cite{Roux99,Kirkwood34,Latimer39,Kornyshev97,Hildebrandt04,Rizzo06,Abrashkin07,Gong09,Guo13,Zhou14}.
Although implicit-solvent models can be orders of magnitude faster
than explicit-solvent molecular-dynamics (MD) simulations, most
popular continuum theories ignore numerous potentially important
effects, including solvent molecules' finite size and specific
molecular interactions such as hydrogen bonding (the AGBNP2 model,
which addresses the latter, is a notable
exception~\cite{Gallicchio09}).

One of the Poisson model's most perplexing and long-standing
shortcomings is the difficulty of extending it to model charge-sign
asymmetric solvation: for example, given two monatomic ions of equal
radius, one of $+q$ charge and the other of $-q$, the negative charge
experiences stronger interactions with the solvent (more negative
solvation free
energy)~\cite{Latimer39,Rashin85,Ashbaugh00,LyndenBell01,Rajamani04,Grossfield05,Mukhopadhyay12,Bardhan12_asymmetry}.
However, standard Poisson models are charge-sign \textit{symmetric};
that is, they predict the same solvation free energy for $\pm q$.  The
need to include asymmetric effects is difficult to exaggerate,
particularly in biological contexts.  Consider that the protein avidin
binds its ligand biotin with a binding free energy of approximately
$-$20~kcal/mol, one of the most favorable in biology~\cite{Green75};
solvent-exposed $+1e$ and $-1e$ charges can experience as much as
40~kcal/mol difference in their solvation free
energies~\cite{Bardhan12_asymmetry}.  Dominant factors in charge-sign
asymmetric response include the liquid-vapor interface
potential~\cite{Ashbaugh00,Rajamani04} and the fact that water
hydrogens can approach a negative solute charge closer than water
oxygens can approach a positive
one~\cite{Latimer39,Rashin85,Mukhopadhyay12,Bardhan12_asymmetry}.
Spherical solutes with central charges provide a useful data set for
developing an understanding of size- and charge-sign dependent
hydration, including the characterization of interface potentials,
solvent packing, and dielectric
saturation~\cite{Rashin85,Hummer96_netchargecorrection,Ashbaugh00,Rajamani04,Grossfield05,Fedorov07}.
These analyses and the continuum macroscopic-dielectric framework
suggest that improvements require a more detailed, accurate
representation of the solvent dipole field
$\mathbf{P}(\mathbf{r})$~\cite{Warshel76,Azuara08}, or, equivalently,
the solvent charge density ~\mbox{$\rho_{\mathrm{induced}}(\mathbf{r})
  = \nabla \cdot \mathbf{P}(\mathbf{r})$}.  Because $\mathbf{P}$ and
$\rho_{\mathrm{induced}}$ do \textit{not} respond linearly to the
solute charge distribution~\cite{Alper90}, particularly in the first
solvent shell~\cite{Purisima09}, many groups have developed solvent
models in which the solvent potential obeys a nonlinear partial
differential equation
(PDE)~\cite{Kornyshev98,Sandberg02_first_LD,Jha08,Gong08,Hu12_Wei_nonlinear_Poisson_BJ}.
Unfortunately, most of these models are still charge-sign symmetric.

However, in 1939 Latimer et al. proposed an approach to increase or
decrease an ion's radius based on the charge~\cite{Latimer39}, and
recent developments in high-performance computing and explicit-solvent
MD free-energy calculations provide important new data to extend this
approach.  Mobley et al. constructed a challenging test set and
conducted extensive MD simulations on charge-sign
asymmetry~\cite{Mobley08_asymmetry}, enabling important new
developments in modeling
asymmetry~\cite{Purisima09,Corbeil10,Mukhopadhyay14} that extend
Latimer's work to Generalized-Born (GB) models of complex solutes.  GB
theory was a natural setting for these developments because Latimer's
work and GB theory share the conceptual picture of an effective atomic
radius.  These early studies provided an important insight: a buried
charge still affects the electric field at the boundary, so merely
parameterizing charge-dependent radii cannot (indeed, should not)
provide a satisfactory explanation.  The accuracy of asymmetric GB
models suggests that a simple Poisson-based model exists, but finding
one has proven to be surprisingly difficult.

%Purisima et al. extended
%Latimer's work to complex solutes using boundary-integral equations
%(BIE), adjusting the radii of solvent-exposed atoms using the local
%induced-charge density
%$\sigma(\mathbf{r})$~\cite{Purisima09,Corbeil10}.  The Onufriev
%group's model interpolates between the limiting cases for zero
%interatomic distances (self-energies) and large
%distances~\cite{Mukhopadhyay14}.

In this paper, we propose a simple Poisson continuum model that
includes charge-sign asymmetry and show that it is remarkably accurate
even without parameterization on an atom-by-atom basis.  The key
feature of our theory is a \textit{nonlinear boundary condition}
(NLBC) for the normal displacement field; in contrast, the
displacement boundary condition for the standard (symmetric) Poisson
theory is linear~\cite{Jackson_classical_electrodynamics}.
Importantly, even though our proposed displacement boundary condition
is nonlinear, the electrostatic potential in the solvent and solute
volumes still satisfy \textit{linear} Poisson/Laplace equations.  Two
phenomena motivated us to propose a nonlinear boundary condition
instead of a nonlinear governing equation.  First, numerous results
illustrate that the solute reaction potential obeys nearly linear
response even though the solvent charge distribution does
not~\cite{Lin11_Pettitt,Boda09,Purisima09,Corbeil10,Mukhopadhyay14,Bardhan12_asymmetry}.
For example, the new asymmetric Generalized-Born (GB) models use the
charge distribution only to modify the Born radii, with the overall
energy still computed using superposition (independent sum of
individual charge
responses)~\cite{Purisima09,Corbeil10,Mukhopadhyay14}.  Furthermore,
we found in our previous work that the solute reaction potential is
essentially a \textit{piecewise-linear} function of charge
~\cite{Bardhan12_asymmetry}, i.e. the proportionality coefficient
depends on whether one is charging an ion from zero to $+q$ or to
$-q$.  In fact, we began this work seeking primarily to reproduce this
curiously simple nonlinearity.

The second phenomenon motivating our NLBC approach is the fact that
the solute reaction potential is a harmonic field---that is, it
satisfies the Laplace equation.  This property is useful for numerical
computations~\cite{Chern03,Holst12} and also provides a path to
improve models via boundary-integral methods~\cite{Bardhan12_review}:
harmonicity means that regardless of the solvent model of interest,
there exists \textit{some} surface charge density that reproduces the
reaction potential inside.  For a given solvent model, the surface
charge density might satisfy a nonlinear boundary-integral equation,
but the very fact that such a density always exists suggests that one
might improve continuum models by adding nonlinear terms to widely
used BIE formalisms~\cite{Rizzo67,Shaw85,Juffer91,Bardhan09_disc}. 

%The following section presents the standard linear model from a
%boundary-integral equation (BIE) perspective, its extension to a
%nonlinear boundary condition (NLBC), and the numerical implementation.
%Section~\ref{sec:results} presents computational results, and the
%discussion in Section~\ref{sec:discussion} addresses the model's
%strengths and weaknesses, as well as implications for implicit-solvent
%models.

\section{Continuum Model and Extension to Nonlinear Boundary Conditions}
We first present the standard (charge-sign symmetric) Poisson
electrostatic model and then describe the difference between it and
our proposed NLBC model.  In both theories, the molecular solute is
treated as a macroscopic linear dielectric continuum obeying the
Poisson equation \mbox{$\nabla^2 \varphi_1(\mathbf{r}) =
  -\frac{\rho(\mathbf{r})}{\epsilon_1}$} where $\mathbf{r}$ is a point
in space, $\varphi_1(\mathbf{r})$ is the potential in the solute,
$\epsilon_1$ is the relative permittivity, and the molecular charge
distribution $\rho(\mathbf{r})$ is a set of $N_q$ point charges, i.e.
\mbox{$\rho(\mathbf{r}) = \sum_{i=1}^{N_q} q_i
  \delta(\mathbf{r}-\mathbf{r}_i)$}.  The solute and solvent are
separated by the interface $\Gamma$, and the solvent exterior is a
linear dielectric with permittivity \mbox{$\epsilon_2 \gg
  \epsilon_1$}, so the electric potential obeys \mbox{$\nabla^2
  \varphi_2(\mathbf{r}) = 0$}; note that modeling realistic biological
solutions requires inclusion of screening effects due to mobile ions
using e.g. some form of the Poisson--Boltzmann equation for
$\varphi_2(\mathbf{r})$~\cite{Sharp90,Tomasi94}.  From macroscopic
dielectric theory and Gauss's law, we obtain the standard Maxwell
boundary conditions for $\mathbf{r}_\Gamma \in \Gamma$
\begin{align}
  \varphi_1(\mathbf{r}_\Gamma) &= \varphi_2(\mathbf{r}_\Gamma),\\
  \epsilon_1 \frac{\partial \varphi_1}{\partial n}(\mathbf{r}_\Gamma) &= \epsilon_2 \frac{\partial \varphi_2}{\partial n}(\mathbf{r}_\Gamma),\label{eq:standard-Maxwell-bc}
\end{align}
where $\frac{\partial}{\partial n}$ denotes the normal derivative (the
normal at $\mathbf{r}_\Gamma$ is defined pointing outward into
solvent). Assuming that $\varphi_2(\mathbf{r})$ decays sufficiently
quickly as $|\mathbf{r}|\rightarrow\infty$, this mixed-dielectric
Poisson problem is well posed and the unknown potential $\varphi_1$
can be rewritten as a linear boundary-integral equation for an unknown
surface charge distribution on $\Gamma$.  In particular, the
apparent-surface charge (ASC)
model~\cite{Shaw85,Altman05_2,Bardhan09_disc} (also known as the
polarizable continuum model~\cite{Miertus81,Tomasi94}) can be
interpreted as finding an equivalent surface charge
$\sigma(\mathbf{r})$ in a homogeneous medium with permittivity
$\epsilon_1$ everywhere.  In this equivalent problem, the
analogous boundary condition to Eq.~\ref{eq:standard-Maxwell-bc}
is simpler due to homogeneity, but adds a term for the surface charge:
\begin{equation}
\frac{\sigma(\mathbf{r}_\Gamma)}{\epsilon_1} = \frac{\partial
    \hat{\varphi}_1}{\partial n}(\mathbf{r}_\Gamma) - \frac{\partial
    \hat{\varphi}_2}{\partial n}(\mathbf{r}_\Gamma),\label{eq:ecf-bc}
\end{equation}
and we use $\hat{\varphi}_i = \varphi_i$ to emphasize our use of an
equivalent problem.  Defining \mbox{$G(\mathbf{r};\mathbf{r}') =
  \frac{1}{4 \pi || \mathbf{r} - \mathbf{r}'||}$}, one obtains
\begin{align}
  \left(I  + \hat{\epsilon} \left(-\frac{1}{2} I+ K\right)\right)\sigma &= -\hat{\epsilon}\sum_i^{N_q} q_i \frac{\partial G}{\partial n}
\end{align}
where $\hat{\epsilon} = (\epsilon_2-\epsilon_1)/\epsilon_2$ and $K$ is
the normal electric field operator~\cite{Bardhan09_disc}.  The
reaction potential in the solute is then
\mbox{$\varphi^{REAC}(\mathbf{r}) = \frac{1}{\epsilon_1} \int_\Gamma
  G(\mathbf{r}; \mathbf{r}') \sigma(\mathbf{r}') dA'$}, and
\mbox{$\varphi_{1}(\mathbf{r}) = \varphi^{REAC}(\mathbf{r})
  +\varphi^{Coulomb}(\mathbf{r})$}, with the latter term representing
the Coulomb potential due to $\rho(\mathbf{r})$.

The standard Maxwell displacement boundary condition
Eq.~\ref{eq:standard-Maxwell-bc} is obtained using Gauss's law in
integral form and the fact that the divergence of the polarization
field $\mathbf{P}(\mathbf{r})$ represents a volume charge density.
However, near the solute--solvent boundary, the assumption that
$\mathbf{P}(\mathbf{r})$ is pointwise proportional to the local
electric field breaks down due to water structure at the interface;
that is, it is no longer necessarily true that
\mbox{$\mathbf{P}(\mathbf{r}) = (\epsilon(\mathbf{r})-1)
  \mathbf{E}(\mathbf{r})$}.

To model nonlinear solvent response at the boundary, we propose to
replace the linear boundary condition,
Eq.~\ref{eq:standard-Maxwell-bc}, with the phenomenological nonlinear
boundary condition
\begin{equation}
f(E_n) \frac{\partial \varphi_1}{\partial n}(\mathbf{r}_\Gamma) = \left(1+f(E_n)\right)\frac{\partial \varphi_2}{\partial n}(\mathbf{r}_\Gamma)\label{eq:charge-layer-nonlinear-BC}
  \end{equation}
where $E_n$ is the electric field just inside $\Gamma$,
i.e. \mbox{$E_n = -\sum_i q_i \frac{\partial G}{\partial n} - K
  \sigma$}, and
\begin{align}
f(E_n) & = \frac{\epsilon_1}{\epsilon_2-\epsilon_1} - h(E_n);\\
h(E_n) & = \alpha \tanh(\beta E_n - \gamma) + \mu.\label{eq:tanh}
  \end{align}
with $\alpha$, $\beta$, and $\gamma$ representing model parameters and
\mbox{$\mu = -\alpha
  \tanh(-\gamma)$}. The specification of $\mu$ ensures that $h(E_n=0)
= 0$, so that in the limit of weak electric fields, such as induced at
the surface by a deeply buried charge, the boundary condition reduces
to the familiar Poisson model.  The NLBC leads to
the modified, \textit{nonlinear} BIE
\begin{equation}
\left(I + \hat{\epsilon}\left(-\frac{1}{2}I + K\right) + h(E_n)\right) \sigma = -\hat{\epsilon}\sum_{i} q_i \frac{\partial G}{\partial n},
  \end{equation}
with the nonlinearity arising in the dependence of $h$ on $E_n$ (see
the Supporting Information for details on the numerical
implementation).

One challenge in developing more accurate solvent models is the fact
that nonlinear response~\cite{Sharp90_2} generally requires a charging
process~\cite{Zhou94}, i.e. the expression \mbox{$\Delta G^{solv,es} = \frac{1}{2}
  q^T \varphi^{REAC} = \frac{1}{2}q^T L q$} no longer holds ($L$
denotes the reaction-potential
operator~\cite{Roux99,Bardhan12_review}). However, our previous work
showed the remarkable fact that the solute reaction potential is
\textit{piecewise} linear, with the breakpoint at
$q=0$~\cite{Bardhan12_asymmetry}, so that $\varphi^{REAC} = L_+ q$ for
$q > 0$ and $\varphi^{REAC} = L_- q$ for $q < 0$, with $L_+ \neq L_-$.
The proposed NLBC in Eqs.~\ref{eq:charge-layer-nonlinear-BC}
and~\ref{eq:tanh} immediately explains this curious phenomenon:
consider the limit $\beta \rightarrow \infty$, so that $\tanh$ is
constant everywhere, but discontinuous at $q=0$.  The Debye charging
process~\cite{Kirkwood34_2} scales all charges uniformly, i.e.
\mbox{$\hat{\rho}(r; \lambda) = \lambda \rho(r)$}, so the Coulomb
field \mbox{$\frac{\partial \varphi^{Coul}}{\partial n}$} has the same
sign for all finite $\lambda$. The Coulomb-field approximation (CFA)
shows that the reaction field is nearly proportional to the direct
Coulomb field, but slightly smaller in
magnitude~\cite{Kharkats76,Bardhan08_BIBEE,Bardhan12_review}, so for
finite $\lambda$, at almost all $\mathbf{r}_\Gamma$, the total field
$E_n(\mathbf{r}_\Gamma; \lambda)$ has the same sign as
$E_n(\mathbf{r}_\Gamma; \lambda = 1)$.  This implies that almost
everywhere on the surface, the $\tanh$ boundary condition takes its
limiting ($\lambda \rightarrow 1$) value for any finite $\lambda$,
which means that the boundary condition is essentially linear:
\mbox{$(1 + g(r))\sigma(r) = \frac{\partial \hat\varphi_2}{\partial
    n}-\frac{\partial \hat\varphi_1}{\partial n}$}.  With this
justification, in this work we compute solvation free energies as
\mbox{$\Delta G^{solv,es}= \frac{1}{2} q^T \varphi^{REAC}$}.  Note
that a more precise definition of the charging free energy would be
piecewise \textit{affine}, because the charging free energy also
includes a linear term that results from the liquid-vapor interface
potential~\cite{Harder08,Kathmann11,Bardhan12_asymmetry}; as noted
above, however, in the present work its influence is approximated via
the offset parameter $\gamma$.

\section{Results and Discussion}\label{sec:results}

We parameterized the NLBC model using the Mobley et al. MD free-energy
calculations, who studied asymmetry using fictitious bracelet and rod
molecules~\cite{Mobley08_asymmetry} constructed from AMBER
$\mathrm{C}_{\alpha}$ atoms with $\mathrm{R}_{\mathrm{min}}/2 = 1.908$~\AA.  We
obtained optimal results with $\alpha = 0.5$, $\beta = -60$, $\gamma =
-0.5$, and a continuum-model $\mathrm{C}_\alpha$ radius of 1.75~\AA~(a
scale factor of approximately 0.92). Note that in this first
exploration of the NLBC, we have parameterized against the overall
solvation free energies computed by Mobley et al. rather than the more
correct charging free energy.

Figure~\ref{fig:standard-ions-model} plots NLBC and MD free-energy
calculations for ion charging free energies; the MD charging
simulations used in our previous work~\cite{Bardhan12_asymmetry} (see
Supporting Information) and CHARMM Lennard-Jones parameters.  We
remind the reader that no additional parameters were fit in obtaining
these NLBC results, i.e. ion radii were assigned
$\mathrm{R}_{\mathrm{ion}} = 0.92 \mathrm{R}_{\mathrm{min}}/2$.  For
additional data, ions were charged to both $+1e$ and $-1e$, regardless
of the charge on the real ion, and the NLBC accurately predicts these
charging free energies as well.  The largest deviations occur for
radii less than 1.4~\AA, where discrete packing effects and actual
dielectric saturation are likely.
\begin{figure}[ht!]
  \centering \resizebox{3.0in}{!}{\includegraphics{vary-ion-radius}}
  \caption{Asymmetric polarization free energies for a monovalent
    central charge in a sphere, as a function of sphere radius.  The
    labeled symbols denote results from MD free-energy calculations
    charging CHARMM monatomic ions from zero to $+1e$ or $-1e$, with
    $\mathrm{R}_{\mathrm{ion}} = 0.92 \mathrm{R}_{\mathrm{min}}/2$.  The dashed black curve
    in the middle is the (charge-sign symmetric) Born polarization
    free energy.}\protect\label{fig:standard-ions-model} \end{figure}

\begin{table}
  \centering
\begin{tabular}{l|cc|cc}
  Problem & \multicolumn{2}{c}{Solvation errors} &  \multicolumn{2}{c}{Asymmetry errors}\\\hline
  & RMSE &  Max.  & RMS & Max. \\
  Rods & 5.57 & 9.63 & 0.88 & 1.49  \\
  Bracelets (opposing) & 2.88 & 6.10 & 2.04 & 3.08  \\
  Bracelets (distributed) & 2.20 & 2.72 & 0.29 & 0.59 \\
  Bracelets (dipole) & 2.67 & 3.52 & 0.85 & 1.09 
\end{tabular}
\caption{Comparison of NLBC model to MD free-energy calculations of
  Mobley et al.~\cite{Mobley08_asymmetry} for rod and bracelet
  molecule test set.  All energies are in kcal/mol.  See Supporting
  Information for detailed results.}\protect\label{table:Mobley}
\end{table}
Calculations for the Mobley test set are summarized in
Table~\ref{table:Mobley}; SI Figures~1--8 plot the NLBC and Mobley MD
solvation free energies and asymmetry energies, and include
illustrations of the test problems.  The rod molecules are composed of
5 or 6 atoms along a line, with one atom possessing $+1e$ charge, one
with $-1e$, and the rest neutral.  The asymmetry errors in
Table~\ref{table:Mobley} represent the difference in solvation
energies when reversing the charged atoms' signs.  The bracelet
molecules are regular polygons with between 3 and 8 sides; atoms are
at the vertices (1.4~\AA~apart). Bracelets were simulated with three
charge distributions: the ``opposing'' case had a $+1e$ charge
neutralized by two $-0.5e$ charges positioned symmetrically on the
opposite side.  The ``distributed'' case has one $+1e$ charge and a
neutralizing $-1e$ distributed equally on all the other atoms; the
``dipole'' case is similar to ``opposing'' but fixes the dipole
moment~\cite{Mobley08_asymmetry}.  Solvent charge-densities from the
MD calculations~\cite{Mobley08_asymmetry} suggest that solvent packing
may be responsible for size-dependent deviations; parameterizing radii
for actual atoms should significantly reduce these errors.

%We reiterate that these calculations used the same atomic radius
%$\mathrm{C}_\alpha = 1.75$\AA~for all test cases.  Small systematic
%deviations in the distributed and dipole cases can be corrected by
%increasing the $\mathrm{C}_\alpha$ radius without changing the NLBC
%parameters, but this increases errors in the opposing cases and the
%rod molecules.

%For the Mobley rod molecules (SI Figures~1 and~2), the NLBC
%predictions have an RMSD error of 5.57~kcal/mol (max. 9.63~kcal/mol)
%and the RMSD error in the solvation-energy asymmetry is 0.88~kcal/mol
%(max. 1.49~kcal/mol).  The Mobley bracelet molecules consist of
%regular polygons with 1.4~\AA~separation between adjacent atoms (SI
%Figures 3--8), and the MD calculations addressed three types of charge
%distributions.  The ``opposing'' charge distributions are the most
%challenging, ranging in solvation free energy from $-$8.19 to
%$-$59.1~kcal/mol (SI Figures~3 and~4).  For these cases the NLBC
%exhibits an RMSD error of 2.88~kcal/mol (max. 6.10~kcal/mol) for
%solvation free energies, with errors exhibiting a similar size
%dependence as in the rod problems. The RMSD error for the asymmetry is
%2.04~kcal/mol (max. 3.08~kcal/mol). In the ``distributed'' charge
%cases, the neutralizing charge is spread equally across the other
%atoms (SI Figures~5 and~6), and the ``dipole'' cases fix the overall
%molecular dipole moment (SI Figures~7 and~8).  
%
%The distributed cases
%have an RMSD error of 2.20~kcal/mol in solvation free energy
%(max. 2.72~kcal/mol), and the dipole cases have an RMSD of
%2.67~kcal/mol (max. 3.52~kcal/mol).  The asymmetries in these cases
%are reproduced very well: for the distributed cases the RMSD error is
%0.29~kcal/mol (max. 0.59~kcal/mol); for the dipole cases, the RMSD
%error is 0.85~kcal/mol (max. 1.09~kcal/mol).  

To test the model on real but nonspherical molecules, we compared NLBC
and MD charging free energies for isolated titratable amino acids in
both protonated and unprotonated states (See Supporting Information
for details on structure preparation). Parameters were from the CHARMM
force field~\cite{Brooks83} when available, with other protonation
states defined so that the protonated and unprotonated states had the
same number of atoms. The MD free-energy-perturbation (FEP)
calculations used the same protocol as the
ions~\cite{Bardhan12_asymmetry}, holding the solute rigid so that
$\epsilon_1 = 1$ unambiguously~\cite{Roux99}.  The deviations between
our MD results and the MD calculations of Nina et al.~\cite{Nina97}
are small compared to the energies of interest, and likely due to our
use of (i) periodic boundary conditions, (ii) a larger solvent box
(1959 waters vs. 150), and (iii) slightly different backbone angles.

As in the ion and Mobley examples, the NLBC radii were defined by the
scaling $\mathrm{R} = 0.92 \mathrm{R}_{\mathrm{min}}/2$. The results
in Figure~\ref{fig:residues} illustrate that the NLBC model correctly
captures solvation free energies in both charge states, despite the
fact that radii were not adjusted individually or even for the atomic
charges.  In contrast, standard Poisson model results computed using
the Nina et al.~\cite{Nina97} or PARSE~\cite{Sitkoff94} radii exhibit
larger deviations, particularly for arginine, aspartic acid, cysteine,
glutamic acid, and tyrosine.  These data suggest that the differences
between symmetric and asymmetric electrostatic models are robust with
respect to radii (the PARSE calculations are merely suggestive because
these calculations used the CHARMM charges; for consistent comparison
to experiment, one should use PARSE charges with PARSE radii).
\begin{figure}[ht!]
  \centering \resizebox{3.0in}{!}{\includegraphics{residues}}
  \caption{Comparison of NLBC model to explicit-solvent MD FEP
    calculations for titratable residues with neutral blocking
    groups. MD results from Nina et al.~\cite{Nina97} are shown where
    available; standard continuum model results are shown for Nina et
    al. radii and PARSE radii~\cite{Sitkoff94} (using CHARMM
    charges).}\protect\label{fig:residues}
\end{figure}

\section{Conclusion}\label{sec:conclusion}

We have proposed a Poisson-based theory that models
charge-sign-dependent asymmetries in electrostatic solvation free
energies using a nonlinear boundary condition (NLBC), while still
using linear continuum theory in the solute and solvent volumes.  The
NLBC model accurately reproduces MD free-energy results for monatomic
ions, the Mobley et al. bracelet and rod problems, and titratable
residues, even though we have used charge-independent radii that were
fixed by a single scaling factor applied to MD radii.  Furthermore,
the NLBC reduces smoothly to the standard Poisson model as the
parameter $\alpha$ approaches zero.  Finally, our boundary-element
method implementation for non-trivial molecules demonstrates that the
new model is easily implemented in numerical Poisson and
Poisson--Boltzmann solvers.

Our introduction of a modified boundary condition to account for
solvation-shell response follows a long history in continuum
mechanics, where phenomenological techniques find applications in many
areas of science and engineering to capture a particular physical
behavior in continuum theory rather than modeling or deriving it from
first
principles~\cite{MizziBarberEmersonReeseStefanov2007,Brenner2011,Bardhan11_Knepley}. Non-equilibrium
micro-scale gas flows offer a well-developed example: velocity-slip
and temperature-jump boundary conditions are simplified
phenomenological approaches to represent both non-equilibrium and
gas-surface interaction effects occurring near solid walls. Such
boundary conditions were first suggested in the 19th century by
Maxwell~\cite{Maxwell1878} and von
Smoluchowski~\cite{Smoluchowski1898}, respectively.  More recent
examples include the partitioning of minerals at phase boundaries in
geophysics~\cite{LHeureux96}, tumor growth~\cite{Macklin05}, the
deformation of biological membranes~\cite{Fan03}, and thin electric
double layers in electro-osmotic flow~\cite{Yossifon07}.

Much as Beglov and Roux showed that solvent response approaches the
linear Poisson model in the limit as the solvent molecule approaches
zero size~\cite{Beglov96}, our model emphasizes that the nonlinear
response is generally localized in the first solvent shell.
Conceptually, the nonlinear boundary condition penalizes negative
surface charge because the larger water oxygen cannot approach a
solute charge as closely as the water hydrogens can.  From a
boundary-integral point of view, this has the same effect as adjusting
the atomic radii, an approach pioneered by Latimer et
al.~\cite{Latimer39}, and extended recently to GB
models~\cite{Purisima09,Corbeil10,Mukhopadhyay12,Mukhopadhyay14}.
Purisima's work is particularly relevant due to their use of
surface-charge boundary-integral approach, adjusting GB radii using
$\sigma(\mathbf{r})$~\cite{Purisima09,Corbeil10}.  Our work differs
substantially from these approaches because we have included asymmetry
directly in the underlying Poisson model.

The present theory can be extended in several important ways.  First,
the proposed NLBC model has only three parameters whose particular
dependencies on solvent model have not yet been established
theoretically.  Second, it seems straightforward to include ionic
screening via the Poisson--Boltzmann equation.  Third, the proposed
NLBC depends exclusively on the normal electric field; improved models
might include local curvature or higher-order moments of the
potential.  Importantly, the latter could distinguish between
small-magnitude charges near the surface, and larger charges further
away~\cite{Mukhopadhyay12}.  Fourth, water's length-scale-dependent
dielectric behavior might be included using nonlocal
electrostatics~\cite{Hildebrandt04,Fedorov07,Bardhan11_pka,Bardhan11_DAC,Bardhan13_nonlocal_review}. The
new model also does not necessarily capture specific hydrogen-bonding
effects like AGBNP2 does~\cite{Gallicchio09}, which motivates future
work comparing the two approaches.

%Future work will explore the possibility of modeling dielectric
%saturation via additional $\tanh$ functions centered at the
%appropriate electric field strengths.

\iftoggle{fulltitlepage}
{
\section*{Acknowledgments}
The authors thank David Mobley for sharing detailed calculation
results, Jed Brown and David Green for valuable discussions, and Matt
Reuter for a critical reading of the manuscript. MGK was partially
supported by the U.S. Department of Energy, Office of Science,
Advanced Scientific Computing Research, under Contract
DE-AC02-06CH11357, and also NSF Grant OCI-1147680. JPB has been
supported in part by the National Institute of General Medical
Sciences (NIGMS) of the National Institutes of Health (NIH) under
award number R21GM102642. The content is solely the responsibility of
the authors, and does not necessarily represent the official views of
the National Institutes of Health.

\section*{Supporting Information}
Figures comparing NLBC and MD calculations for the Mobley test
set~\cite{Mobley08_asymmetry}.  The source code (MATLAB) and surface
discretizations for running the nonlinear boundary-condition
calculations, data files, parameters, and scripts for preparing and
running the MD calculations of titratable residues, as well as source
code to generate the figures, are freely and publicly available online
at \url{https://bitbucket.org/jbardhan/si-nlbc}.
}{
}

\bibliographystyle{unsrt}
\bibliography{implicit-review}

\end{document}


\section*{Extra Results}

Points for reader
\begin{enumerate}
 \item Atomic radii have not been fit to optimize.  We used RminOver2
  for the ion results and the Mobley examples (for the latter, the
  atoms used $\sigma = 3.39967$~\AA, so $Rmin/2 = (2^{1/6} \sigma)/2 =
  1.908$~\AA, and for the 5A sphere we used the central charge case
  (obtaining 5.06 and 5.11 for +1 and -1) where asymmetry does not
  play a role.
\item We're talking about the electrostatic charging free energy
  only. The bracelet and rod cases are total hydration energies, and
  these cases are challenging with significant non-bulk and
  non-electrostatic components.
  \end{enumerate}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Points about Figure~\ref{fig:vary-ion-radius}
\begin{enumerate}
\item Details on our FEP calculations may be found in earlier
  work~\cite{Bardhan12_asymmetry} or in the supplemental information.
\item The results in the Figure do not include the contribution from
  the liquid-vapor interface potential (static
  potential)~\cite{Bardhan12_asymmetry}.  
\item Much work on continuum models has focused on adjusting ion
  radii~\cite{Latimer39,Rashin85}.  Here, we have not suggested
  ion radii for the continuum model, but merely plotted the continuum
  model predictions together with the results from MD/FEP, using
  CHARMM Rmin/2 parameters as possible radii.
\item The present NLBC model is a simple ansatz and cannot account for
  real dielectric saturation; for example, for polyvalent ions such as
  calcium, the first solvation shell is maximally polarized and then
  solvent response proceeds in the second shell.
\item Deviations may be seen even in lithium, for example.
\end{enumerate}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[ht!]
  \centering \resizebox{6.0in}{!}{\includegraphics{../figures/compare-bjm-sphere}}
  \caption{Charging free energies for a charge in a 5.1~\AA~sphere.
    The symbols denote results from MD free-energy
    calculations~\cite{Bardhan12_asymmetry}.  Both the MD and the NLBC model
    energies include the contribution from the liquid-vapor interface
    potential (static potential)~
    \cite{Bardhan12_asymmetry}.}\protect\label{fig:compare-bjm-sphere} \end{figure}

Points about Figure~\ref{fig:compare-bjm-sphere}
\begin{enumerate}
  \item The sphere was constructed by creating a Cartesian grid
    (spacing 0.5~\AA) of uncharged atoms with sodium's Lennard-Jones
    parameters, so that no atom was $\geq 4$~\AA~from the origin.
  \item Other details, including about the MD calculations, can be
    found in~\cite{Bardhan12_asymmetry}
  \item To determine a radius for this nearly spherical molecule, we
    began with the free energy for the +1$e$ central charge, removed
    the static potential (10.4~kcal/mol/$e$) contribution, and then
    inverted the Born equation for the radius (5.1~\AA).  Nearly
    identical results are obtained using the -1$e$ central charge
    (5.06~\AA).
  \item The central charge experiences very little steric asymmetry
    (approx. 1 kcal/mol), cf. Fig.~\ref{fig:vary-ion-radius} where the
    steric asymmetry is very small for even a 4~\AA~sphere.  
  \item Our model therefore correctly captures the dependence of
    solvation free energy on the degree of burial.
  \end{enumerate}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[ht!]
  \centering \resizebox{6.0in}{!}{\includegraphics{../figures/Mobley-fig9-b}}
  \caption{}\protect\label{fig:Mobley-fig9-b} \end{figure}

Points about Figures~\ref{fig:Mobley-fig9-a} and~\ref{fig:Mobley-fig9-b}
\begin{enumerate}
\item We generated structures according to the descriptions of Mobley et al.~\cite{Mobley08}.
\item Also following their work, our continuum calculations used Bondi
  radii for all of the atoms (1.7~\AA).
\item Similar to the ring results, our energies usually are too
  negative, meaning that the radius could be increased. However,
  fitting a radius for fictitious molecules seems like a waste of
  time.  We will, however, need to parameterize for biomolecules.
\item The absolute solvation free energies vary significantly, ranging
  between -15~kcal/mol and -110~kcal/mol.  However, the magnitude of
  the asymmetry changes very little, and our model reproduces this
  very well, even though NLBC consistently underestimates the
  asymmetry by approximately 2~kcal/mol.
\item Notably, the semi-buried ion pair (the charges at positions 2
  and 3 in the 5-member rod) experiences a reversal in asymmetry,
  which our model reproduces.  This may be fortuitous, but a main
  point is that all of these charges are solvent exposed at least on
  the sides.
\item Although the physical interpretation of changes from one case to
  another is not necessarily meaningful, the correlation between the
  results seems good; the notable exceptions are the greatest charge
  separations for the 6-member case, where solvent discreteness
  effects can play a larger role.
\item we might be concerned about the overall magnitude being
  off. note that if you change alpha to zero, we get the same
  magnitude, so parameterizing around solvation free energies should
  restore us to a model with both solvation and asymmetry right
\end{enumerate}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[ht!]
  \centering \resizebox{6.0in}{!}{\includegraphics{../figures/Mobley-fig5-c}}
  \caption{}\protect\label{fig:Mobley-fig5-c} \end{figure}
\begin{figure}[ht!]
  \centering \resizebox{6.0in}{!}{\includegraphics{../figures/Mobley-fig5-d}}
  \caption{}\protect\label{fig:Mobley-fig5-d} \end{figure}
\begin{figure}[ht!]
  \centering \resizebox{6.0in}{!}{\includegraphics{../figures/Mobley-fig4-a}}
  \caption{}\protect\label{fig:Mobley-fig4-a} \end{figure}
\begin{figure}[ht!]
  \centering \resizebox{6.0in}{!}{\includegraphics{../figures/Mobley-fig4-b}}
  \caption{}\protect\label{fig:Mobley-fig4-b} \end{figure}

Points about Figures~\ref{fig:Mobley-fig5-c}-\ref{fig:Mobley-fig4-b}:
\begin{enumerate}
\item Again structures and radius taken from Mobley~\cite{Mobley08}
\item Our predictions consistently underestimate the asymmetry by 2-3~kcal/mol
\end{enumerate}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[ht!]
  \centering \resizebox{6.0in}{!}{\includegraphics{../figures/Mobley-fig5-a}}
  \caption{}\protect\label{fig:Mobley-fig5-a} \end{figure}
\begin{figure}[ht!]
  \centering \resizebox{6.0in}{!}{\includegraphics{../figures/Mobley-fig5-b}}
  \caption{}\protect\label{fig:Mobley-fig5-b} \end{figure}

Points about Figures~\ref{fig:Mobley-fig5-a} and \ref{fig:Mobley-fig5-b}:
\begin{enumerate}
\item Again structures and radius taken from Mobley~\cite{Mobley08}
\item 
\end{enumerate}


\section*{References to add}
\begin{enumerate}
\item \cite{Fedorov07} Kornyshev/Fedorov review of the nonlocal models
\item \cite{Harder08} Roux interface potential asymmetry
\item \cite{Kornyshev97} Kornyshev+Sutmann nonlocal dielectric saturation
\item idea: if you had a double layer source at the surface couldn't you also generate the static potential effect?  what would that total boundary condition look like?  would it be phiInt = phiExt + doubleLayerSource?
  \item New Dill and Roux papers on interface potential
  \item \cite{Paluch00,Cerutti07,Harder08,Kathmann11} liquid-vapor potential interface
  \item \cite{Schutz01} Schutz Warshel review on implicit structural relaxation
  \item \cite{Gunner00} backbone dipoles positive potential in proteins
  \end{enumerate}

\section{Basic Derivation}
\begin{align}
  \varphi_1 &= \varphi_2\\
  \frac{\partial \varphi_1}{\partial n} & =  \frac{\partial \varphi_2}{\partial n}\\
  \nabla^2 \varphi_1 &= - \frac{\rho(r)}{\epsilon_1} \\
  \rho(r) & =  \sum_{i=1}^{N_c} q_i \delta(r-r_i)\\
  \nabla^2 \varphi_2 &= 0\\
  \frac{\sigma(r)}{\epsilon_1} &= \frac{\partial \hat{\varphi}_2}{\partial n} - \frac{\partial \hat{\varphi}_1}{\partial n} \\
  G(r;r') & =  \frac{1}{4 \pi \epsilon_1 || r - r'||}\\
  \hat{\varphi}_1(r) &= \sum_{i=1}^{N_c} q_i G(r; r') + \int_\Gamma \sigma(r') G(r; r') dA'\\
  \frac{\partial \hat{\varphi}}{\partial n}(r_\Gamma) &= \sum_{i=1}^{N_c} q_i \frac{\partial}{\partial n} G(r; r') + \frac{\alpha}{2} \sigma + \dashint_\Gamma \sigma(r') \frac{\partial}{\partial n} G(r; r') dA'\\
  & = \sum q_i \frac{\partial G}{\partial n} + \left(\frac{\alpha}{2\epsilon_1}I + K\right) \sigma
  \end{align}
I never remember whether the limit gives $+1/2 I$ or $-1/2 I$ so I
leave it as $\alpha$ for now and we'll double check later.  Assuming
that we're taking the limit from the inside, i.e. in the $\varphi_1$
region, then we proceed by writing $\varphi_2$ in terms of $\varphi_1$
\begin{align}
  \frac{\partial \hat{\varphi}_2}{\partial n}& = \frac{\epsilon_1}{\epsilon_2} \frac{\partial \hat{\varphi}_1}{\partial n}\\
  \frac{\sigma}{\epsilon_1} &= \frac{\epsilon_1}{\epsilon_2}\frac{\partial \hat{\varphi}_1}{\partial n} - \frac{\partial \hat{\varphi}_1}{\partial n}\\
  & = \frac{\epsilon_1 - \epsilon_2}{\epsilon_2}\frac{\partial \hat{\varphi}_1}{\partial n}\\
  \frac{\epsilon_2}{\epsilon_1(\epsilon_1-\epsilon_2)}\sigma &=\frac{\partial \hat{\varphi}_1}{\partial n}\\
  & = \sum q_i \frac{\partial G}{\partial n} + \left(\frac{\alpha}{2\epsilon_1}I+K\right) \sigma\\
  \left(\frac{\epsilon_2}{\epsilon_1(\epsilon_1-\epsilon_2)}-\frac{\alpha}{2\epsilon_1}-K\right)\sigma &= \sum q_i \frac{\partial G}{\partial n}\\
  \left(\frac{\alpha}{2 \epsilon_1} - \frac{\epsilon_2}{\epsilon_1(\epsilon_1-\epsilon_2)} + K\right)\sigma &= -\sum q_i \frac{\partial G}{\partial n}
\end{align}
so if $\alpha = +1/2$
\begin{align}
\frac{\alpha}{2 \epsilon_1} - \frac{\epsilon_2}{\epsilon_1(\epsilon_1-\epsilon_2)} & =
 \frac{(\epsilon_1-\epsilon_2) - 2 \epsilon_2}{2 \epsilon_1 (\epsilon_1 - \epsilon_2)}\\
 & = \frac{\epsilon_1 - 3 \epsilon_2}{2 \epsilon_1(\epsilon_1 - \epsilon_2)}
  \end{align}
and if $\alpha = -1/2$
\begin{align}
\frac{\alpha}{2 \epsilon_1} - \frac{\epsilon_2}{\epsilon_1(\epsilon_1-\epsilon_2)} & =
\frac{-(\epsilon_1-\epsilon_2) - 2 \epsilon_2}{2 \epsilon_1 (\epsilon_1-\epsilon_2)} \\
& = \frac{-(\epsilon_1 + \epsilon_2)}{2 \epsilon_1(\epsilon_1 - \epsilon_2)}\\
& = \frac{\epsilon_1 + \epsilon_2}{2 \epsilon_1(\epsilon_2-\epsilon_1)}\\
& = \frac{\bar{\epsilon}}{\epsilon_1 \Delta \epsilon}
  \end{align}
where $\bar{\epsilon} = \frac{1}{2}(\epsilon_1+\epsilon_2)$ and
$\Delta \epsilon = \epsilon_2 - \epsilon_1$ (i.e. the mean and the
difference crossing over in the normal direction).  This reveals an
unfortunate sign error in the Altman EMBC paper, but at least the
signs are right in the PRE paper with Dirk, and in our other BIBEE
papers.  So $\alpha = -1/2$.



Problems with existing theories
\begin{enumerate}
\item Charges near the surface that are nonetheless buried? they
  experience asymmetric solvation (see our paper with Pavel) but many
  existing models only get asymmetry right for surface charge
  (i.e. buried charges ``affect the radius'' of the surface atoms, but
  this is not included in many models).  One of the few models that
  DOES get this right is Onufriev's CHA-GB.
\item Models need to include the static potential.
\end{enumerate}


Discussion points
\item Relatively easy to implement in large scale solvers (BC only
  adjusted, not PDE; Picard converges quickly, given Coulomb E field
  starting guess, so we do not need to do Newton,), both volumetric
  and boundary-integral

\item improved models could account for?  (i.e. not drive diagonal
  entry to change sign: intuitively unphysical)

\item divalent ions show different breakpoints. can this be included? should this be in discussion?

\item Model works by total normal electric field at interface, but simpler (not-self-cons) version can be obtained by focusing on bare-Coulomb field from solute charges

\item Nonlinear boundary condition can be included in both PDE and BIE approaches easily

\item simple model that invokes nonlinear boundary conditions but enforces discontinuous dielectric
\item \cite{Alper90} Alper and Levy field strength dependence of dielectric saturation in water

\item Models should take care to account for the static potential
  separately because it's actually linear wrt charge.

\item Also extensible to PCM models

draft.bbl
\begin{thebibliography}{10}

\bibitem{Roux99}
B.~Roux and T.~Simonson.
\newblock Implicit solvent models.
\newblock {\em Biophys. Chem.}, 78:1--20, 1999.

\bibitem{Sharp90}
K.~A. Sharp and B.~Honig.
\newblock Electrostatic interactions in macromolecules: Theory and
  applications.
\newblock {\em Annu. Rev. Biophys. Bio.}, 19:301--332, 1990.

\bibitem{Tomasi94}
J.~Tomasi and M.~Persico.
\newblock Molecular interactions in solution: {An} overview of methods based on
  continuous descriptions of the solvent.
\newblock {\em Chem. Rev.}, 94:2027--2094, 1994.

\bibitem{Kirkwood34}
J.~G. Kirkwood.
\newblock Theory of solutions of molecules containing widely separated charges
  with special application to zwitterions.
\newblock {\em J. Chem. Phys.}, 2:351, 1934.

\bibitem{Latimer39}
W.~M. Latimer, K.~S. Pitzer, and C.~M. Slansky.
\newblock The free energy of hydration of gaseous ions, and the absolute
  potential of the normal calomel electrode.
\newblock {\em J. Chem. Phys.}, 7:108--112, 1939.

\bibitem{Kornyshev97}
A.~A. Kornyshev and G.~Sutmann.
\newblock Nonlocal dielectric saturation in liquid water.
\newblock {\em Phys. Rev. Lett.}, 79:3435--3438, 1997.

\bibitem{Hildebrandt04}
A.~Hildebrandt, R.~Blossey, S.~Rjasanow, O.~Kohlbacher, and H.-P. Lenhof.
\newblock Novel formulation of nonlocal electrostatics.
\newblock {\em Phys. Rev. Lett.}, 93:108104, 2004.

\bibitem{Rizzo06}
R.~C. Rizzo, T.~Aynechi, D.~A. Case, and I.~D. Kuntz.
\newblock Estimation of absolute free energies of hydration using continuum
  methods: accuracy of partial charge models and optimization of nonpolar
  contributions.
\newblock {\em J. Chem. Theory Comput.}, 2:128--139, 2006.

\bibitem{Abrashkin07}
A.~Abrashkin, D.~Andelman, and H.~Orland.
\newblock Dipolar {Poisson}--{Boltzmann} equation: ions and dipoles close to
  charge interfaces.
\newblock {\em Phys. Rev. Lett.}, 99:077801, 2007.

\bibitem{Gong09}
H.~Gong and K.~F. Freed.
\newblock {Langevin--Debye} model for nonlinear electrostatic screening of
  solvated ions.
\newblock {\em Phys. Rev. Lett.}, 102(057603), 2009.

\bibitem{Rashin85}
A.~A. Rashin and B.~Honig.
\newblock Reevaluation of the {Born} model of ion hydration.
\newblock {\em J. Phys. Chem.}, 89:5588--5593, 1985.

\bibitem{Ashbaugh00}
Henry~S. Ashbaugh.
\newblock Convergence of molecular and macroscopic continuum descriptions of
  ion hydration.
\newblock {\em The Journal of Physical Chemistry B}, 104(31):7235--7238, 2000.

\bibitem{LyndenBell01}
R.~M. {Lynden-Bell}, J.~C. Rasaiah, and J.~P. Noworyta.
\newblock Using simulation to study solvation in water.
\newblock {\em Pure Appl. Chem.}, 73:1721--1731, 2001.

\bibitem{Rajamani04}
S.~Rajamani, T.~Ghosh, and S.~Garde.
\newblock Size dependent ion hydration, its asymmetry, and convergence to
  macroscopic behavior.
\newblock {\em J. Chem. Phys.}, 120:4457, 2004.

\bibitem{Grossfield05}
A.~Grossfield.
\newblock Dependence of ion hydration on the sign of the ion's charge.
\newblock {\em J. Chem. Phys.}, 122:024506, 2005.

\bibitem{Mukhopadhyay12}
A.~Mukhopadhyay, A.~T. Fenley, I.~S. Tolokh, and A.~V. Onufriev.
\newblock Charge hydration asymmetry: the basic principle and how to use it to
  test and improve water models.
\newblock {\em J. Phys. Chem. B}, 116:9776–9783, 2012.

\bibitem{Bardhan12_asymmetry}
J.~P. Bardhan, P.~Jungwirth, and L.~Makowski.
\newblock Affine-response model of molecular solvation of ions: Accurate
  predictions of asymmetric charging free energies.
\newblock {\em J. Chem. Phys.}, 137:124101, 2012.

\bibitem{Green75}
N.~M. Green.
\newblock Avidin.
\newblock {\em Adv. Protein Chem.}, 29:85--133, 1975.

\bibitem{Hummer96_netchargecorrection}
G.~Hummer, L.~R. Pratt, and A.~E. Garc\'{i}a.
\newblock Free energy of ionic hydration.
\newblock {\em J. Phys. Chem.}, 100:1206--1215, 1996.

\bibitem{Fedorov07}
M.~V. Fedorov and A.~A. Kornyshev.
\newblock Unravelling the solvent response to neutral and charged solutes.
\newblock {\em Molecular Physics}, 105:1--16, 2007.

\bibitem{Warshel76}
A.~Warshel and M.~Levitt.
\newblock Theoretical studies of enzymic reactions: Dielectric, electrostatic
  and steric stabilization of the carbonium ion in the reaction of lysozyme.
\newblock {\em J. Mol. Biol.}, 103:227--249, 1976.

\bibitem{Azuara08}
C.~Azuara, H.~Orland, M.~Bon, P.~Koehl, and M.~Delarue.
\newblock Incorporating dipolar solvents with variable density in
  {Poisson}--{Boltzmann} electrostatics.
\newblock {\em Biophys. J.}, 95:5587--5605, 2008.

\bibitem{Alper90}
H.~E. Alper and R.~M. Levy.
\newblock Field strength dependence of dielectric saturation in liquid water.
\newblock {\em J. Phys. Chem.}, 94:8401--8403, 1990.

\bibitem{Purisima09}
E.~O. Purisima and T.~Sulea.
\newblock Restoring charge asymmetry in continuum electrostatic calculations of
  hydration free energies.
\newblock {\em J. Phys. Chem. B}, 113:8206--8209, 2009.

\bibitem{Kornyshev98}
A.~A. Kornyshev and G.~Sutmann.
\newblock Nonlocal nonlinear static dielectric response of polar liquids.
\newblock {\em J. Electroanalytical Chemistry}, 450:143--156, 1998.

\bibitem{Sandberg02_first_LD}
L.~Sandberg and O.~Edholm.
\newblock Nonlinear response effects in continuum models of the hydration of
  ions.
\newblock {\em J. Chem. Phys.}, 116:2936--2944, 2002.

\bibitem{Jha08}
A.~K. Jha and K.~F. Freed.
\newblock Solvation effect on conformations of {1,2:dimethoxyethane}:
  Charge-dependent nonlinear response in implicit solvent models.
\newblock {\em J. Chem. Phys.}, 128(034501), 2008.

\bibitem{Gong08}
H.~Gong, G.~Hocky, and K.~F. Freed.
\newblock Influence of nonlinear electrostatics on transfer energies between
  liquid phases: Charge burial is far less expensive than {Born} model.
\newblock {\em Proc. Natl. Acad. Sci. USA}, 105:11146--11151, 2008.

\bibitem{Hu12_Wei_nonlinear_Poisson_BJ}
L.~Hu and G.-W. Wei.
\newblock Nonlinear {Poisson} equation for heterogeneous media.
\newblock {\em Biophys. J.}, 103:758--766, 2012.

\bibitem{Mobley08_asymmetry}
D.~L. Mobley, A.~E. {Barber II}, C.~J. Fennell, and K.~A. Dill.
\newblock Charge asymmetries in hydration of polar solutes.
\newblock {\em J. Phys. Chem. B}, 112:2405--2414, 2008.

\bibitem{Corbeil10}
C.~R. Corbeil, T.~Sulea, and E.~O. Purisima.
\newblock Rapid prediction of solvation free energy. 2. {The} first-shell
  hydration {(FiSH)} continuum model.
\newblock {\em J. Chem. Theory Comput.}, 6:1622--1637, 2010.

\bibitem{Mukhopadhyay14}
A.~Mukhopadhyay, B.~H. Aguilar, I.~S. Tolokh, and A.~V. Onufriev.
\newblock Introducing charge hydration asymmetry into the {Generalized Born}
  model.
\newblock {\em J. Chem. Theory Comput.}, 10:1788--1794, 2014.

\bibitem{Jackson_classical_electrodynamics}
J.~D. Jackson.
\newblock {\em Classical Electrodynamics}.
\newblock Wiley, 3$^{rd}$ edition, 1998.

\bibitem{Lin11_Pettitt}
B.~Lin, K.-Y. Wong, C.~Hu, H.~Kokubo, and B.~M. Pettitt.
\newblock Fast calculations of electrostatic solvation free energy from
  reconstructed solvent density using proximal radial distribution functions.
\newblock {\em J. Phys. Chem. Lett.}, 2:1626--1632, 2011.

\bibitem{Boda09}
D.~Boda, M.~Valisk\'{o}, D.~Henderson, D.~Gillespie, B.~Eisenberg, and M.~K.
  Gilson.
\newblock Ions and inhibitors in the binding site of {HIV} protease:
  {Comparison} of {Monte Carlo} simulations and the linearized
  {Poisson}--{Boltzmann} theory.
\newblock {\em Biophys. J.}, 96:1293--1306, 2009.

\bibitem{Chern03}
I.-L. Chern, J.-G. Liu, and W.-C. Wang.
\newblock Accurate evaluation of electrostatics for macromolecules in solution.
\newblock {\em Methods Appl. Anal.}, 10:309--328, 2003.

\bibitem{Holst12}
M.~Holst, J.~A. {McCammon}, Z.~Yu, Y.~C. Zhou, and Y.~Zhu.
\newblock Adaptive finite element modeling techniques for the
  {Poisson}--{Boltzmann} equation.
\newblock {\em Commun. Comput. Phys.}, 11:179--214, 2012.

\bibitem{Bardhan12_review}
J.~P. Bardhan.
\newblock Biomolecular electrostatics---{I} want your solvation (model).
\newblock {\em Computational Science and Discovery}, 5:013001, 2012.

\bibitem{Rizzo67}
F.~J. Rizzo.
\newblock An integral equation approach to boundary value problems of classical
  electrostatics.
\newblock {\em Quarterly of Applied Mathematics}, 25:83--95, 1967.

\bibitem{Shaw85}
P.~B. Shaw.
\newblock Theory of the {Poisson Green's}-function for discontinuous dielectric
  media with an application to protein biophysics.
\newblock {\em Phys. Rev. A}, 32(4):2476--2487, 1985.

\bibitem{Juffer91}
A.~H. Juffer, E.~F.~F. Botta, B.~A.~M. van Keulen, A.~van~der Ploeg, and
  H.~J.~C. Berendsen.
\newblock The electric potential of a macromolecule in a solvent: A fundamental
  approach.
\newblock {\em J. Comput. Phys.}, 97(1):144--171, 1991.

\bibitem{Bardhan09_disc}
J.~P. Bardhan.
\newblock Numerical solution of boundary-integral equations for molecular
  electrostatics.
\newblock {\em J. Chem. Phys.}, 130:094102, 2009.

\bibitem{Altman05_2}
M.~D. Altman, J.~P. Bardhan, J.~K. White, and B.~Tidor.
\newblock An efficient and accurate surface formulation for biomolecule
  electrostatics in non-ionic solution.
\newblock In {\em Engineering in Medicine and Biology Conference (EMBC)}, 2005.

\bibitem{Miertus81}
S.~Miertus, E.~Scrocco, and J.~Tomasi.
\newblock Electrostatic interactions of a solute with a continuum -- a direct
  utilization of \textit{ab initio} molecular potentials for the prevision of
  solvent effects.
\newblock {\em Chem. Phys.}, 55(1):117--129, 1981.

\bibitem{Sharp90_2}
K.~A. Sharp and B.~Honig.
\newblock Calculating total electrostatic energies with the nonlinear
  {Poisson}--{Boltzmann} equation.
\newblock {\em J. Phys. Chem.}, 94(19):7684--7692, 1990.

\bibitem{Kirkwood34_2}
J.~G. Kirkwood.
\newblock On the theory of strong electrolyte solutions.
\newblock {\em J. Chem. Phys.}, 2:767--772, 1934.

\bibitem{Kharkats76}
{Yu.}~I. Kharkats, A.~A. Kornyshev, and M.~A. Vorotyntsev.
\newblock Electrostatic models in the theory of solutions.
\newblock {\em J. Chem. Soc. Faraday Trans.2}, 72:361--371, 1976.

\bibitem{Bardhan08_BIBEE}
J.~P. Bardhan.
\newblock Interpreting the {Coulomb}-field approximation for
  {Generalized}-{Born} electrostatics using boundary-integral equation theory.
\newblock {\em J. Chem. Phys.}, 129(144105), 2008.

\bibitem{Harder08}
E.~Harder and B.~Roux.
\newblock On the origin of the electrostatic potential difference at the
  liquid-vapor interface.
\newblock {\em J. Chem. Phys.}, 129:234706, 2008.

\bibitem{Kathmann11}
S.~M. Kathmann, I-F.~W. Kuo, C.~J. Mundy, and G.~K. Schenter.
\newblock Understanding the surface potential of water.
\newblock {\em J. Phys. Chem. B}, 115:4369--4377, 2011.

\bibitem{Brooks83}
B.~R. Brooks, R.~E. Bruccoleri, B.~D. Olafson, D.~J. States, S.~Swaminathan,
  and M.~Karplus.
\newblock {CHARMM:} {A} program for macromolecular energy, minimization, and
  dynamics calculations.
\newblock {\em J. Comput. Chem.}, 4:187--217, 1983.

\bibitem{Nina97}
M.~Nina, D.~Beglov, and B.~Roux.
\newblock Atomic radii for continuum electrostatics calculations based on
  molecular dynamics free energy simulations.
\newblock {\em J. Phys. Chem. B}, 101:5239--5248, 1997.

\bibitem{Sitkoff94}
D.~Sitkoff, K.~A. Sharp, and B.~Honig.
\newblock Accurate calculation of hydration free energies using macroscopic
  solvent models.
\newblock {\em J. Phys. Chem. B}, 98:1978--1988, 1994.

\bibitem{MizziBarberEmersonReeseStefanov2007}
S.~Mizzi, R.W. Barber, D.R.Emerson, J.M. Reese, and S.K. Stefanov.
\newblock A phenomenological and extended continuum approach for modelling
  non-equilibrium flows.
\newblock {\em Continuum Mech. Thermodyn.}, 19:273--283, 2007.

\bibitem{Brenner2011}
Howard Brenner.
\newblock Beyond the no-slip boundary condition.
\newblock {\em Physical Review E}, 84(4), 2011.

\bibitem{Bardhan11_Knepley}
J.~P. Bardhan and M.~G. Knepley.
\newblock Mathematical analysis of the boundary-integral based electrostatics
  estimation approximation for molecular solvation: Exact results for spherical
  inclusions.
\newblock {\em J. Chem. Phys.}, 135:124107, 2011.

\bibitem{Maxwell1878}
James~Clerk Maxwell.
\newblock On stresses in rarefied gases arising from inequalities of
  temperature.
\newblock {\em Proceedings of the Royal Society of London}, 27:304--308, 1878.

\bibitem{Smoluchowski1898}
Marian von Smolan~Smoluchowski.
\newblock {\"Uber} w\"armeleitung in verd\"unnten gasen.
\newblock {\em Annalen der Physik}, 300(1):101--130, 1898.

\bibitem{LHeureux96}
I.~{L'Heureux} and A.~D. Fowler.
\newblock Dynamical model of oscillatory zoning in plagioclase with nonlinear
  partition relation.
\newblock {\em Geophysical Research Letters}, 23:17--20, 1996.

\bibitem{Macklin05}
P.~Macklin and J.~Lowengrub.
\newblock Evolving interfaces via gradients of geometry-dependent interior
  {Poisson} problems: application to tumor growth.
\newblock {\em J. Comput. Phys.}, 203:191--220, 2005.

\bibitem{Fan03}
T.-H. Fan and A.~G. Fedorov.
\newblock Electrohydrodynamics and surface force analysis in {AFM} imaging of a
  charged, deformable biological membrane in a dilute electrolyte solution.
\newblock {\em Langmuir}, 19:10930--10939, 2003.

\bibitem{Yossifon07}
G.~Yossifon, I.~Frankel, and T.~Miloh.
\newblock Symmetry breaking in induced-charge electro-osmosis over polarizable
  spheroids.
\newblock {\em Physics of Fluids}, 19:068105, 2007.

\bibitem{Beglov96}
D.~Beglov and B.~Roux.
\newblock Solvation of complex molecules in a polar liquid: an integral
  equation theory.
\newblock {\em Journal of Chemical Physics}, 104(21):8678--8689, 1996.

\bibitem{Bardhan11_pka}
J.~P. Bardhan.
\newblock Nonlocal continuum electrostatic theory predicts surprisingly small
  energetic penalties for charge burial in proteins.
\newblock {\em J. Chem. Phys.}, 135:104113, 2011.

\bibitem{Bardhan11_DAC}
J.~P. Bardhan and A.~Hildebrandt.
\newblock A fast solver for nonlocal electrostatic theory in biomolecular
  science and engineering.
\newblock In {\em IEEE/ACM Design Automation Conference (DAC)}, 2011.

\bibitem{Bardhan13_nonlocal_review}
J.~P. Bardhan.
\newblock Gradient models in molecular biophysics: progress, challenges,
  opportunities.
\newblock {\em Journal of Mechanical Behavior of Materials}, 22:169--184, 2013.

\bibitem{Jayaram89}
B.~Jayaram, R.~Fine, K.~Sharp, and B.~Honig.
\newblock Free energy calculations of ion hydration: an analysis of the {Born}
  model in terms of microscopic simulations.
\newblock {\em J. Phys. Chem.}, 93:4320--4327, 1989.

\end{thebibliography}
