# Golden aspect ratio for ion transport simulation in nanopores

Subin Sahu  Center for Nanoscale Science and Technology, National Institute of
Standards and Technology, Gaithersburg, MD 20899  Maryland NanoCenter,
University of Maryland, College Park, MD 20742  Michael Zwolak<EMAIL_ADDRESS>
Center for Nanoscale Science and Technology, National Institute of Standards
and Technology, Gaithersburg, MD 20899

###### Abstract

Access resistance indicates how well current carriers from a bulk medium can
converge to a pore or opening, and is an important concept in nanofluidic
devices and in cell physiology. In simplified scenarios, when the bulk
dimensions are infinite in all directions, it depends only on the resistivity
and pore radius. These conditions are not valid in all-atom molecular dynamics
(MD) simulations of transport, due to the computational cost of large
simulation cells, and can even break down in micro- and nano-scale systems due
to strong confinement. Here, we examine a scaling theory for the access
resistance that predicts a special simulation cell aspect ratio – the golden
aspect ratio – where finite size effects are eliminated. Using both continuum
and all-atom simulations, we demonstrate that this golden aspect ratio exists
and that it takes on a universal value in linear response and moderate
concentrations. Outside of linear response, it gains an apparent dependence on
characteristics of the transport scenario (concentration, voltages, etc.) for
small simulation cells, but this dependence vanishes at larger length scales.
These results will enable the use of all-atom molecular dynamics simulations
to study contextual properties of access resistance – its dependence on
protein and molecular-scale fluctuations, the presence of charges, and other
functional groups – and yield the opportunity to quantitatively compare
computed and measured resistances.

##  I INTRODUCTION

Ion transport through narrow constrictions in biological membranes permits the
regulation of concentrations, known as ion homeostasis, that is vital for
physiological functions of cells Hille (2001); Kandel  _et al._ (2000); Doyle
_et al._ (1998). Moreover, ion transport through porous inorganic membranes is
of interest for technological applications – such as DNA sequencing Zwolak and
Di Ventra (2008); Heerema and Dekker (2016) – and industrial use Baker (2000);
Mulder (2012). The recent progress in the fabrication of atomically thin
membranes such as graphene Garaj  _et al._ (2010); Merchant  _et al._ (2010);
Schneider  _et al._ (2010), MoS2 Heiranian  _et al._ (2015) and hexagonal
boron nitride Walker  _et al._ (2017) opens new avenues in the field of ion
transport: These membranes are excellent candidates for molecular and ionic
sieves Joshi  _et al._ (2014); Abraham  _et al._ (2017); Esfandiar  _et al._
(2017) for desalination Cohen-Tanugi and Grossman (2012); Surwade  _et al._
(2015) and gas separation Jiang  _et al._ (2009); Kim  _et al._ (2013). Their
atomic thickness provides advantages for bio-sensing, such as sequencing
Clarke  _et al._ (2009); Schneider  _et al._ (2010); Merchant  _et al._ (2010)
and protein folding Si and Aksimentiev (2017). Synthetic pores – especially
ones with atomic thickness – also provide a testing ground for understanding
biological ion channels and creating biomimetic pores. In particular, ion
currents through pores of controllable size can probe dehydration Zwolak  _et
al._ (2009, 2010) and directly quantify its effect on selectivity Sahu  _et
al._ (2017); Sahu and Zwolak (2017).

Figure 1:  Contour plots of potential from continuum (PNP) simulations of ion
transport through a nanopore. (a) Potential drop $V_{0}$ across the simulation
cell with a pore of radius $a$ in a membrane of thickness $h_{p}$. The
equipotential surfaces (thick, shaded red and blue lines shown on top of a
heat map) are elliptical near the pore – i.e., an access-like region. The
finite size of the simulation cell curtails the access region, forcing the
equipotential surfaces to transition into nearly flat profiles closer to the
electrodes – i.e., a bulk-like region. The top half shows the rotational
elliptical coordinates, $\xi$ and $\eta$, we employ in modeling the access
resistance. (b) Equipotential surfaces in simulation cells of increasing $L$
and $H$, with $H>L$, show the finite-size scaling of access-resistance. (c)
Potential map in simulation cells with $L>H$ ($L=32$ nm and $H=16$ nm). The
equipotential surfaces are quasi-elliptical almost up to the end of the cell,
but with the surfaces somewhat vertically distorted when they approach the
electrode surfaces on the top and bottom. (d) As the cell size increases at a
constant aspect ratio $\alpha=H/L$, more of the equipotential surfaces become
ellipsoidal, i.e., access-like and the total resistance converges to its value
in an infinite bulk.

The access resistance – part of the series resistance in patch clamp
measurements – is the resistance for ions to converge from the bulk
electrolyte to the mouth of the pore Hall (1975). It sets the upper limit of
current flow through ion channels Hille (1968); Läuger (1976) and can become
the dominating resistance at low salt concentration Alcaraz  _et al._ (2017).
For atomically thin pores with radii sufficiently larger than the solvation
shell of ions, the access resistance will be the only significant resistance,
even at high salt concentration. It thus becomes an important component to
understand and quantify for sensing and sequencing, which require high
precision measurement and analysis of the ion current. Moreover, numerous
efforts, especially using water-soluble polymers, have sought to separate pore
and access contributions to the resistance in order to characterize different
aspects of biological channels Vodyanoy and Bezrukov (1992); Bezrukov and
Vodyanoy (1993); Bezrukov  _et al._ (1996).

Ideally, all-atom molecular dynamics (MD) simulations should be employed to
simulate ion channels Maffeo  _et al._ (2012), as only these simulations
capture contextual aspects of the pores, such as molecular-scale fluctuations
in pore sizes/geometries and local charges Sahu and Zwolak (2018). For
instance, edge fluctuations and the noncircular nature of a pore in graphene,
together with van der Waals interactions and dehydration, make the pore radius
hard to define Sahu and Zwolak (2018). However, MD simulations are
computationally intensive. In fact, it is often not possible to reach
biologically relevant timescales using all-atom MD Eisenberg (2010). Together
with the long-range nature of convergence, this makes it difficult to reach
the required simulation sizes to quantify the access resistance Yoo and
Aksimentiev (2015). Here, we examine a scaling analysis to extract the access
and pore resistances Sahu and Zwolak (2018), demonstrating that there is
indeed a “golden aspect ratio” (different than the golden ratio,
$(1+\sqrt{5})/2$) that removes finite size effects. We use continuum
simulations in order to scan the necessary parameter space and demonstrate
that the golden aspect ratio indeed removes finite size effects in MD
simulation as well.

For an infinitely large, balanced (i.e., infinite in all directions), and
homogeneous system, the access resistance depends only on the resistivity of
the medium, $\gamma$, and the pore radius, $a$. Hall’s expression,

$R_{\mathrm{Hall}}=\frac{\gamma}{4a},$ (1)

gives the form under these idealized conditions Hall (1975). In fact, the
access resistance, Eq. 1, was derived by Maxwell well over a century ago in
the context of the electrical current through an orifice Maxwell (1881).
Access resistance occurs in heat flow, mass diffusion, and other related
scenarios, such as electrical (e.g., disc electrodes) and thermal contacts.
Hence, several authors have derived the same expression for these various
physical systems, see, e.g., Gray and Mathews Gray and Mathews (1895), Brown
and Escombe Brown and Escombe (1900), Gröber Gröber (1921), Holm Holm (1958),
and Newmann Newman (1966), and given various names for it (access,
convergence, contact, a component of the series resistance, etc.). In the case
of electronic transport, the resistance in the ballistic regime – the Sharvin
resistance Sharvin (1965) – crosses over to Maxwell’s expression when the pore
is larger than the electron mean free path Wexler (1966); Nikolić and Allen
(1999). Although our focus is on ion transport, the general findings should be
applicable to other transport scenarios as well.

Moreover, the fact that the access resistance varies as $1/a$ rather than
$1/a^{2}$, as it does for the pore resistance in the diffusive regime, has an
interesting consequence: When the pore resistance is negligible, the current
through one large pore is less than that through several smaller pores of the
same total area if the pores do not interfere with each other. Nature uses
this effect to maximize the gas exchange rate between the atmosphere and
stomata in leaves Brown and Escombe (1900), and one can envision using the
same effect for maximizing permeation through atomically thin membranes, such
as those through graphene and MoS2, which naturally have a small pore
resistance (above the dehydration limit Sahu  _et al._ (2017); Sahu and Zwolak
(2017)).

As we discuss later, Eq. 1 reflects an idealized situation. Access resistance
can deviate significantly from this form when the pore and the membrane are
charged Aguilella-Arzo  _et al._ (2005); Luchinsky  _et al._ (2009) or when
only one ion species is permeable Läuger (1976); Peskoff and Bers (1988). The
effect of surface charge or concentration gradients on access resistance have
been studied elsewhere Aguilella-Arzo  _et al._ (2005); Peskoff and Bers
(1988). Here, we focus on the effect of the finite simulation cell size.

##  II Models and Methods

In addition to an infinitely large, balanced, and homogeneous bulk, the
typical derivation of Eq. 1 assumes that a hemispherical electrode is at
infinity. This means that at large distances the electric field lines extend
radially outward from the pore/contact. Within the pore, these field lines
have to transition to pointing along the symmetry axis, which we take as the
$z$-axis. These symmetries can be seen from the equipotential surfaces close
to the pore, as shown in Fig. 1. We will retain the ellipsoidal symmetry,
which transitions from circular at the pore to hemispherical at infinity, to
derive the finite-size corrections. These corrections are especially important
in simulations, but can also be relevant to nanofluidic Panday  _et al._
(2016) and microelectromechanical systems (MEMS) Read  _et al._ (2009).

Figure 2:  Variation of the resistance with changes in one of the cell
dimensions. The pore radius is $a=1$ nm, the membrane thickness $h_{p}=1$ nm,
and the resistivity $\gamma=71$ M$\Omega\cdot$nm (taken to match our all-atom
MD value of 1 M KCl from prior results Sahu and Zwolak (2018). (a) The
resistance $R$ versus cross-sectional length $L$ with electrolyte height fixed
at $H=14.5$ nm. For $L\gtrsim H/f$, $R$ is smaller than the actual resistance
for an infinite-size electrolyte, $R_{\infty}$ (horizontal dotted line). The
scaling form still correctly predicts $R_{\infty}$ when fit for small $L$
(solid lines) but to do so it has to display non-monotonic behavior
(extrapolated, dashed lines). The computed data at large $L$, however,
converges to a value below $R_{\infty}$ as given by Eq. 11. We note that
obtaining $f$ from where $R$ crosses $R_{\infty}$ gives about 1.0 for
cylindrical and 1.2 for rectangular cells. This is in close agreement with the
golden aspect ratio we obtain later, with the difference from the exact values
likely due to the strong non-monotonic effects here. (b) The resistance $R$
versus total electrolyte height $H$ when cross-section is fixed at $L=50$ nm.
For $H/f\gtrsim L$, $R$ increases linearly with increasing $H$ as shown by the
fitted solid lines. For $H<L$, $R$ decreases with decreasing $H$ as the
negative correction term $-f^{\prime}\gamma/H$ becomes larger. The factor
$f^{\prime}\approx 0.4$ obtained from fitting is similar to simple estimates
(about 0.3). The standard error of the fits for $f$ and $R_{\infty}$ are less
than 0.2 %. Their fitted values are shown in bold font in the legends.

We have shown previously Sahu and Zwolak (2018) how to modify Eq. 1 for a
finite-size system [while retaining the symmetry of the problem Braunovic  _et
al._ (2006); Newman (1966)]. This setup is easier to solve using rotational
elliptic coordinates, $\xi$ and $\eta$, which relate to cylindrical
coordinates, $z$ and $\rho$, via

$\displaystyle z$ $\displaystyle=a\,\xi\,\eta$ $\displaystyle\rho$
$\displaystyle=a\sqrt{(1+\xi^{2})(1-\eta^{2})}.$ (2)

Laplace’s equation for the electric potential in this coordinate system is
Abramowitz and Stegun (1964)

$\frac{\partial}{\partial\xi}\left[\left(1+\xi^{2}\right)\frac{\partial
V}{\partial\xi}\right]+\frac{\partial}{\partial\eta}\left[\left(1-\eta^{2}\right)\frac{\partial
V}{\partial\eta}\right]=0.$ (3)

The boundary conditions are: (i) a constant potential at the pore mouth ($V=0$
at $\xi=0$), (ii) no radial electric field on the membrane surface ($\partial
V/\partial\eta=0$ at $\eta=0$), and (iii) an ellipsoidal electrode of radius
$l$ ($V=V_{0}$ at $\xi=l/a$). Only condition (iii) is different than that used
for an infinite bulk.

We stress that in a complex system such as ion channel, factors such as the
presence of surface charges and functional groups, concentration gradients,
selective ion transport, and others limit the validity of these boundary
conditions. Even in a simplified system these boundary conditions only
partially hold. In particular, condition (iii) is really a fictitious
electrode placed at the end of the access region. We will use the location as
a fitting parameter, as well as an additional contribution to capture the
effect of the transition region – the region between the access and bulk
regions. Our main focus will be on obtaining the functional dependence of the
resistance on cell dimensions and we will use continuum simulations to both
empirically motivate and validate this dependence.

The solution to Eq. 3 is, see Appendix,

$\frac{V}{V_{0}}=\frac{\tan^{-1}\xi}{\tan^{-1}(l/a)}$ (4)

for the potential. Thus, the access resistance is

$\displaystyle R_{\mathrm{access/2}}$
$\displaystyle=\frac{\gamma\tan^{-1}(l/a)}{2\pi a}$
$\displaystyle=\frac{\gamma}{4a}\left(1-\frac{2a}{\pi
l}+\mathcal{O}\left[\left(\frac{a}{l}\right)^{3}\right]\right)$
$\displaystyle\approx R_{\mathrm{Hall}}\left(1-\frac{2a}{\pi
l}\right)=R_{\mathrm{Hall}}-\frac{\gamma}{2\pi l},$ (5)

where we ignore higher order corrections since they are
$\mathcal{O}\left[(a/l)^{3}\right]$. The notation “access/2” indicates that
this is the access resistance on a single side of the membrane. When the
(simulation) cell size has, e.g., a cross-sectional length of 10 nm, $l$ has
to be less than 5 nm. Hence, the expression shows that the classical form is
off by $\approx$ 13 % for a pore radius of 1 nm. This will get worse for even
moderately larger pore sizes, not to mention the difficulty in applying the
Hall’s form when $a$ is ill-defined.

Figure 3:  Resistance versus $L$ for different aspect ratios $\alpha$, with
$a=1$ nm, $h_{p}=1$ nm, and $\gamma=71$ M$\Omega\cdot$nm. For small $\alpha$,
$dR/dL>0$ and for large $\alpha$, $dR/dL<0$. At special value of $\alpha$ in
between – the “golden aspect ratio”, $dR/dL=0$ and the resistance is constant
(i.e., no finite size effects). Here, the gold line is very close to, but not
quite at, the golden aspect ratio. The solid lines show the fits
$R=\frac{\gamma}{\mathcal{G}L}\left(\alpha-f\right)+R_{\infty}$ for $\alpha>1$
and $R=-\frac{\gamma f^{\prime}}{\alpha L}+R_{\infty}$ for $\alpha<1$, where
$\mathcal{G}=\pi/4$ for cylindrical cells (left panel) and $\mathcal{G}=1$ for
rectangular cells (right panel). The fit parameters are shown in bold. The
dashed lines extrapolate these fits, which match well with the calculations
for large $L$ and yield consistent values for $R_{\infty}$. The performance of
the scaling analysis using small $L$ indicates that the simulation cell sizes
achievable with all-atom MD should be sufficient to find the total resistance
(pore plus access). The error of the fits for the $f$’s and the $R_{\infty}$’s
are about 0.5 % and 0.1 %, respectively (except for $\alpha=0.25$ where the
respective errors are about about 3 % and 0.5 % ).

Scaling Analysis
The simple boundary conditions and ellipsoidal symmetry allow us to derive the
expression, Eq. II, for the access resistance for an idealized finite size
system. However, as mentioned earlier, neither the boundary conditions nor the
symmetry holds exactly in practice. In ion channels, the potential in the pore
and the bulk can be coupled Luchinsky  _et al._ (2009). Nonetheless, such
coupling affects the region near the pore and the boundary conditions can be
taken as approximations. Another important consideration is that simulations
usually have parallel disc electrodes (or homogeneous applied fields) and a
uniform cross-section. Nanopore experiments and patch-clamp measurement of
biological pores have even more complicated arrangements. Furthermore, in
simulations, a rectangular or cylindrical cell are the natural choices.
However, as shown in Fig. 1, the potential still has ellipsoidal symmetry near
the pore and it starts to become flat away from the pore in the vertical
direction. This is also the case in all-atom MD simulations Sahu and Zwolak
(2018). Most importantly, therefore, we have to consider the empirical
observation, see Fig. 1, that there are different electrostatic regions of the
cell. In what follows, we will use this observation to develop a general
scaling form.

In the case where the bulk height $H$ (not including the membrane thickness)
is greater than the cross-sectional length $L$, i.e., $H\gtrsim L$, there are
three regions on each side: (a) an ellipsoidal access-like region extending
from $\xi=0$ to $\xi=l/a$, (b) a flat bulk-like region extending from
approximately $z=l$ to $z=H/2$ (this takes the upper membrane surface to be at
$z=0$), and (c) an intermediate region between the two. When $H\gtrsim L$, the
length $l$ will be some fraction of the cross-sectional length $L$ and the
resistance of the intermediate region should decay as $1/L$. Thus, the total
resistance is

$\displaystyle R$
$\displaystyle=R_{\mathrm{access}}+R_{\mathrm{bulk}}+R_{\mathrm{intermediate}}+R_{\mathrm{pore}}$
(6) $\displaystyle=2\left(R_{\mathrm{Hall}}-\frac{\gamma}{\pi
f_{1}L}\right)+\gamma\left(\frac{H-f_{2}L}{\mathcal{G}L^{2}}\right)+\frac{\gamma
f_{3}}{\mathcal{G}L}+R_{\mathrm{pore}}$

where $\mathcal{G}L^{2}$ is the cross-sectional area of the bulk cell
($\mathcal{G}=1$ for rectangular and $\mathcal{G}=\pi/4$ for cylindrical). We
introduce factors $f_{1}$, $f_{2}$, and $f_{3}$ due to the uncertainty in the
extent/contribution of these regions and the transitory nature of the
boundaries between them. In Eq. 6, the access region ends at $l=f_{1}L/2$ and
we expect $f_{1}$ to be $\mathcal{O}(1)$, i.e., this region encompasses a
substantial fraction of the cell width (when $H\gtrsim L$). After the access
region and some transition region ends in the vertical direction, the normal
bulk-like region begins. It has total height $H-2f_{2}L/2$, where $f_{2}L/2$
is subtracted from each side of the membrane (the membrane thickness is not
included in $H$), and we expect $f_{2}$ to also be $\mathcal{O}(1)$.

When the cell size is infinite in all directions, the resistance reduces to

$R\to R_{\infty}=2R_{\mathrm{Hall}}+R_{\mathrm{pore}}.$ (7)

Thus, from equations 6 and 7, we obtain

$\displaystyle R=$
$\displaystyle\frac{\gamma}{\mathcal{G}}\left(\frac{H}{L^{2}}-\frac{f}{L}\right)+R_{\infty},$
(8)

where

$f=2\mathcal{G}/\pi f_{1}+f_{2}-f_{3}$ (9)

gives a single fitting factor. This relation shows that the finite size – and
confined – correction to the access resistance depends on both height and the
cross-section of the cell. In the context of heat flow Mikic (1967); Cooper
_et al._ (1969), others have shown that, for an infinitely tall cell, access
resistance has functional dependence on $a/L$, which reduces to the classical
form when $L\gg a$.

Equation 8 is the resistance for $H\gtrsim L$ and its development employed
both the derivation of the corrections to Hall’s form and the empirical
observations in Fig. 1. However, Eq. 8 does not necessarily require
$R_{\infty}=2R_{\mathrm{Hall}}+R_{\mathrm{pore}}$, rather only that the
resistance convergences to some $R_{\infty}$. For example, $R_{\infty}$ can
include the rectifying action of a channel Bhattacharya  _et al._ (2011), the
effect of charges Aguilella-Arzo  _et al._ (2005); Luchinsky  _et al._ (2009),
or an ill-defined pore radius Sahu and Zwolak (2018). Conditions such as
these, ones that give an access contribution other than Hall’s form, should
still obey the scaling law so long as the cell is large enough to remove non-
scaling finite size effects. The reason is simple, the convergence is an
algebraically decaying effect, whereas structural fluctuations are completely
local and charges are screened beyond few Debye lengths. This means that after
some distance these features will no longer be felt by the ionic solution. We
will show this explicitly in the case of non-linear response to large
voltages, which is a much more drastic perturbation to the medium than
fluctuations of the pore or the presence of charges.

Figure 4:  Universality of the golden aspect ratio with pore size. (a) The
resistance versus $L$ for cylindrical cells with $\alpha=1.07$, $h_{p}=1$ nm,
$\gamma=71$ M$\Omega\cdot$nm, and various pore radii. The golden aspect ratio
is fairly constant over this range of $a$. The solid lines show the fits
$R=\frac{4\gamma}{\pi L}\left(\alpha-f\right)+R_{\infty}$ and the dashed lines
show the extrapolations – which make good predictions of $R$ at large $L$. The
fit parameters are shown in bold. The errors of the fits for $f$’s and
$R_{\infty}$’s are less than 0.5 % and 0.1 % respectively. As the pore radius
increases, we examine fits for larger $L$, to roughly keep proportionality
between $L$ and $a$ as there are (non-scaling) finite-size effects that will
affect the simulations when the radius starts to become comparable to the cell
cross-sectional length. The range of the fitting can influence the extracted
$f$ (we expect that ultra-precise calculations at really large $L$ will reveal
the exact golden aspect ratio, which will be close to 1.07). (b) $\Delta R$
from fits to $R=\gamma\left(1/2a+h_{\mathrm{eff}}/\pi a^{2}\right)+\Delta R,$
for cylindrical cells of various $\alpha$, with $h_{\mathrm{eff}}$ and $\Delta
R$ as fitting parameters. We get $\Delta R\approx 0$ for the “golden aspect
ratio” $\alpha^{\star}=1.07$ (dashed lines), the fit of which is shown in the
inset. The effective height, $h_{\mathrm{eff}}$, is (5$\pm$1) % larger than
the 1 nm membrane thickness, see Fig. 5. (c) $\Delta R$ in the fits
$R=\gamma/4a+\Delta R$ for half cylinders of different aspect ratios with the
pore mouth set at potential $V=0$. We again get $\Delta R\approx 0$ for
$\alpha^{\star}=1.07$ (dashed lines), whose fit is shown in the inset. The
interpolated lines in the insets are for visual clarity only.

For $H\lesssim L$, we will get a different form than Eq. 8 but the
considerations will be similar. Fig. 1 (c) shows the equipotential surfaces
for this case. The access region now ends at $l=f_{1}^{\prime}H/2$, as it is
not the horizontal boundary that terminates the access region but the
vertical. This occurs at some fraction of the box height,
$f_{1}^{\prime}=\mathcal{O}(1)$, where the prime indicates it is different
than the fraction in the $H\gtrsim L$ case. The normal bulk region is now
almost negligible. Moreover, the cross-sectional area that contributes to the
resistance in the normal bulk region is not $L\times L$, but rather should be
nearly equal to $H\times H$, as only the region above the access region feeds
ions into that region and contributes to the current. This is also reflected
by the fact that the electric field is nearly zero for a radial distance about
$H$ away from the $z$-axis. Hence, at that distance away, the electrolyte does
not contribute to the resistance. We will therefore take $\gamma
f_{2}^{\prime}/\mathcal{G}H$ for the normal bulk contribution (i.e., a height
divided by a relevant cross sectional length that is proportional to
$H/H^{2}$) and $f_{2}^{\prime}$ should be very small. The transition region
will also depend on $1/H$, giving

$\displaystyle R$
$\displaystyle=R_{\mathrm{access}}+R_{\mathrm{bulk}}+R_{\mathrm{intermediate}}+R_{\mathrm{pore}}$
(10) $\displaystyle=2\left(R_{\mathrm{Hall}}-\frac{\gamma}{\pi
f_{1}^{\prime}H}\right)+\frac{\gamma
f_{2}^{\prime}}{\mathcal{G}H}+\frac{\gamma
f_{3}^{\prime}}{\mathcal{G}H}+R_{\mathrm{pore}}.$

Rewriting this equation gives

$\displaystyle R=$ $\displaystyle-\frac{\gamma f^{\prime}}{H}+R_{\infty},$
(11)

with $f^{\prime}=1/\pi
f_{1}^{\prime}-f_{2}^{\prime}/\mathcal{G}-f_{3}^{\prime}/\mathcal{G}$. This
equation entails that if we keep $H$ constant the total resistance will stay
constant as we increase $L$. We estimate $f_{1}^{\prime}\approx 1$ and
$f_{2}^{\prime}\approx f_{3}^{\prime}\approx 0$ giving $f^{\prime}\approx
1/\pi\approx 0.32$. Moreover, the correction to access resistance is negative
as it removes part of the access region to give a lower lower total
resistance.

The difficulty with the applying Eq. 8 and 11 is that $R$ both depends on two
geometric variables, $H$ and $L$, and it transitions from Eq. 8 to Eq. 11 at
an unknown boundary (albeit $H\approx L$). On the one hand, if we hold $L$
constant and increase $H$, then $R$ will increase linearly with $H$ (e.g.,
this would apply both in simulations and experiments where a very narrow
nanofluidic constriction leads up to the membrane/pore). This is essentially
increasing the larger normal bulk region in the cell. On the other hand, if we
keep $H$ constant and increase $L$, the resistance will initially decrease due
to the larger cross-sectional area available for transport in the normal bulk
region. However, the finite-size correction in Eq. 8 becomes negative once
$L\gtrsim H/f$. At this point, $R$ becomes smaller than $R_{\infty}$ and, as
Eq. 11 indicates, should flatten out as $L$ becomes larger. This is seen in
Fig. 2. In other words, a further increase in $L$ does not converge $R$ to
$R_{\infty}$, because the additional cross-sectional area does not contribute
to transport. One can still apply the scaling form, Eq. 8, so long as the fit
is for a range of $L\lesssim H/f$, which Fig. 2(a) shows it extracts the
correct $R_{\infty}$ (found by using the approach we describe below). To use
the scaling form in Eq. 11 one has to set $L$ to a large value and increase
$H$, fitting for values $H\lesssim L$, as seen in Fig. 2(b) (one can not take
$H$ near in magnitude to $L$, however, as there are corrections missing in Eq.
10 as $H$ approaches $L$ from below).

In order for $R$ to converge to $R_{\infty}$ monotonically (in the absence of
nonlinearities), whether from above or from below, both $L$ and $H$ need to
increase simultaneously. An intuitive method to do so is to keep the aspect
ratio, $\alpha=H/L$, constant, which simplifies Eq. 8 and Eq. 11 to

$\displaystyle
R=\begin{cases}\frac{\gamma}{\mathcal{G}L}\left(\alpha-f\right)+R_{\infty}&\text{for
}\alpha\gtrsim 1,\\\ -\frac{\gamma f^{\prime}}{\alpha L}+R_{\infty}&\text{for
}\alpha<1.\end{cases}$ (12)

Not only do these equations give a unified scaling form
$R_{\infty}+\mathcal{O}(1/L)$, i.e., both decay with $L$ (one from above, one
from below), the former also indicates that if we choose a special aspect
ratio $\alpha^{\star}=f$ then the finite size correction is eliminated and $R$
is independent of $L$. In other words, at this “golden aspect ratio”,
$\alpha^{\star}$, $R=R_{\infty}$ for any $L$.

It is likely not possible to find a general expression $f=2\mathcal{G}/\pi
f_{1}+f_{2}-f_{3}$, since the transition from access-like to bulk-like is
complicated and the factors $f_{1}$, $f_{2}$ and $f_{3}$ will in general
depend on geometric details of the pore and the boundary conditions.
Nevertheless, with some reasonable approximations we can give an estimate of
$f$. We expect that $f_{1}\approx f_{2}\approx 1$ (i.e., a boundary between
access and bulk-like regions at a radial distance equal to the cell edge,
$l=L/2$). To estimate $f_{3}$, which depends on $f_{1}$ and $f_{2}$, we assume
that the intermediate region has equal contributions from access- and bulk-
like behavior. This intermediate region can be assumed to have a hemispherical
boundary of radius $f_{1}L/2$ and a flat boundary of cross-section
$\mathcal{G}L^{2}$ at height $f_{2}L/2$ on each side of the pore. Thus the
resistance of the intermediate region is estimated to be

$\frac{\gamma f_{3}}{\mathcal{G}L}=\frac{1}{2}\left(\frac{2\gamma}{\pi
f_{1}L}-\frac{2\gamma}{\pi L}\right)+\frac{1}{2}\left(\frac{2\gamma
f_{2}L/2}{\mathcal{G}L^{2}}\right),$ (13)

giving

$f_{3}=\frac{\mathcal{G}}{f_{1}\pi}-\frac{\mathcal{G}}{\pi}+\frac{f_{2}}{2}.$
(14)

The expression for $f_{3}$ likely overestimates the bulk-like contribution,
while either under- or over-estimating the access contribution. Using it in
Eq. 9, simplifies $f$ to

$\displaystyle f$ $\displaystyle=\frac{\mathcal{G}}{\pi
f_{1}}+\frac{f_{2}}{2}+\frac{\mathcal{G}}{\pi}.$ (15)

$f_{1}$ | 1.00 | 0.80 | 0.60 | 0.40
---|---|---|---|---
$f_{2}$ | 1.00 | 0.80 | 0.60 | 0.40
$\alpha^{\star}$ (rectangular) | 1.14 | 1.12 | 1.15 | 1.31
$\alpha^{\star}$ (cylindrical) | 1.00 | 0.96 | 0.97 | 1.08
Table 1:  Estimates of the golden aspect ratio, $\alpha^{\star}=f$, for
various values of $f_{1}$ and $f_{2}$. We expect $f_{1}$ and $f_{2}$ to be
$\mathcal{O}(1)$.

Table 1 shows the value of $f$ for various values of $f_{2}=f_{1}$. The value
of $f$ is quite insensitive to $f_{1}$ and $f_{2}$ in a reasonable range (due
to the fact that one appears in a denominator and the other in a numerator),
including when $f_{1}$ and $f_{2}$ are varied separately. The equipotential
surfaces in Fig. 1 suggest $f_{1}\approx 2/3$ and $f_{2}\approx 3/4$, and thus
we estimate

$f=\begin{cases}1.2&\text{for rectangular box}\\\ 1.0&\text{for cylindrical
box}\end{cases}$ (16)

for the golden aspect ratio.

We will also examine the sensitivity of $f$ to different conditions, including
whether one uses the Poisson-Nernst-Planck (PNP) equations or just Laplace’s
equation (i.e., Ohm’s law). The latter – the homogeneous medium approximation
– is valid when the applied potential is small and the ion concentration is
large Garaj  _et al._ (2010). Poisson’s equation and the stationary Nernst-
Planck equation express the spatial dependence of the potential $V$ and
current density $\bf{J}$ as

$\nabla^{2}V=\sum_{\nu}\frac{q_{\nu}c_{\nu}}{\epsilon}$ (17)

and

${\bf J_{\nu}}=-q_{\nu}\left(D_{\nu}\nabla c_{\nu}+\mu_{\nu}c_{\nu}\nabla
V\right),$ (18)

where $q_{\nu}$, $c_{\nu}$, $D_{\nu}$, $\mu_{\nu}$ are the charge,
concentration, diffusivity, and mobility of ion species $\nu$. We use a
commercial finite element solver for both cases.

Finally, we expect that the “golden aspect ratio” is applicable to MD because,
at large distance and weak enough variation in fields, MD can be coarse-grain
into a continuum description. We test this applicability by simulating the
ionic current through graphene nanopore of radius $a=1.81$ nm immersed in 1
mol/L of KCl. We perform the all-atom MD simulation using NAMD2 Phillips  _et
al._ (2005), the CHARMM27 Feller and MacKerell (2000) force field, and a rigid
TIP3P Jorgensen  _et al._ (1983) water model. The details of MD simulation are
the same as in our previous work Sahu and Zwolak (2018).

Figure 5:  The resistance versus $a$ for cylindrical cells with $L=32$ nm, the
golden aspect ratio $\alpha^{\star}=1.07$, and $h_{p}=$ 1.0 nm, 2.0 nm, 3.0 nm
and 4.0 nm. The resistance fits into the model,
$R=\gamma\left(1/2a+h_{\mathrm{eff}}/\pi a^{2}\right)+\Delta R$ with
$h_{\mathrm{eff}}$ and $\Delta R$ as fitting parameters (shown in bold in the
legend). The fitted $h_{\mathrm{eff}}$ is slightly larger than the actual
membrane thickness but its relative difference diminishes as the membrane
thickness is increased. This is due to the curvature of the potential at the
pore mouth and not due to numerical grid-size errors (a ten times smaller grid
increase this value to 7 % and thus we expect this height correction to
converge to something of this order.). The standard errors of
$h_{\mathrm{eff}}$ and $\Delta R$ are about 0.01 nm and 0.04 M$\Omega$ in each
case. Figure 6:  PNP solution of the resistance versus $L$ for cylindrical
cells with $h_{p}=1$ nm, $a=1$ nm, and various aspect ratios ($\alpha=$ 1.4,
1.1 and 0.8 from top to bottom). The labels at the top show the applied
potential and the labels at the right show the concentration of KCl. The
golden aspect ratio remains the same as in the homogeneous medium calculation
for the smaller voltages (10 mV and 100 mV) and for both 1 mol/L and 0.1 mol/L
concentrations. For 300 mV applied voltage, the golden aspect ratio seems to
take on a lower value due to the increased concentration of ions near the
membrane, which is more pronounced for low concentrations due to a diminished
ability to screen and generate local fields. The legends show the fits to Eq.
(12) with the fit parameters shown in bold. The error of the fits for $f$’s
and $R_{\infty}$’s are about 0.5 % and 0.1 % respectively.

##  III RESULTS

We first test the scaling form, Eq. 12, for a homogeneous medium by solving
Laplace’s equation for the electric potential. All simulations take a silica
membrane with dielectric constant 2.1. With the rectangular cell, we compute
the resistance using both finite and periodic cells (as is common in all-atom
molecular dynamics). Both boundary conditions give the same results and thus
we show results only for the finite cells for the continuum case (for the all-
atom simulations, we have a periodic rectangular cell). Fig. 3 shows that, for
both cylindrical and rectangular cells, the resistance increases with $L$ for
small aspect ratios and decreases for large aspect ratios as predicted from
Eq. 12. Physically, we can understand this as follows: When the aspect ratio
($H/L$) is small, the access-like region covers most of the height of the
cell. Thus, when increasing $L$, the access contribution increases, as the
finite size correction – a negative correction – in Eq. II is being
eliminated. On the other hand, when the aspect ratio ($H/L$) is large, the
access-like region is localized near the pore and the bulk-like region covers
the remainder of the height of the cell. This bulk-like region gives a large
positive contribution to the resistance, but is decreasing with $L$ (due to
the fact that it decreases with the area). This transition from increasing to
decreasing resistance with $L$ suggests that there should be a constant $R$ at
some special value of the aspect ratio. At this special value, the two
competing effects cancel. Fig. 3 shows that there is an aspect ratio that
shows little variation as $L$ increases. Thus, indeed, the golden aspect ratio
exists and is approximately $\alpha^{\star}=1.07$ for a cylindrical cell and
$\alpha^{\star}=1.2$ for a rectangular cell, in line with the estimate in Eq.
16.

In Fig. 4(a), we examine three pore sizes and use a cylindrical cell with
aspect ratio $\alpha=1.07$ where we get a visually flat profile of $R$ versus
$L$. An inappropriately chosen aspect ratio, however, will have $R$ dependent
on $L$ and disagree with the classical form of the access resistance. Fig.
4(b) shows that, around $\alpha=1.07$ and at finite $L$ (32 nm), the form
$R=\gamma\left(1/2a+h_{p}/\pi a^{2}\right)+\Delta R$, fits the data with
$\Delta R\approx 0$, thus coinciding with the classical form. For $\alpha=2$,
a reasonable value from a computational standpoint, the fit gives $\Delta
R\approx 3$ M$\Omega$, which is an error of $\approx 30$ % in $R_{\infty}$ at
$a=4$ nm.

Notice that we obtain $h_{p}=1.05$ nm from the fit, which is different than
the 1 nm membrane thickness. The expression $\gamma h_{p}/\pi a^{2}$ assumes
the potential surfaces are flat within the full length of the pore. However,
in actual pores, these surfaces are not perfectly flat, especially near the
pore mouth, as is visible in Fig. 1. This difference in curvature will give a
correction to $\gamma h_{p}/\pi a^{2}$ and it seems likely this is responsible
for the slight difference in $h_{p}$. This correction becomes negligible for
thick membranes, as shown in Fig. 5, where the correction comes out to be 1.0
% for a 4 nm membrane (compared to 5 % for the 1 nm membrane).

To test the finite-size effect of the bulk, we can remove this source of
ambiguity completely by examining the resistance on a half cylinder – between
the pore mouth set at potential $V=0$ and an electrode at one end of the
cylinder. Fig. 4(c) shows this resistance versus pore radius. Since the half
cylinder does not have the pore resistance, we fit it to $R=\gamma/4a+\Delta
R$. Once again we get $\Delta R\approx 0$ for the aspect ratio $\alpha=1.07$
and a value substantially different than zero for other $\alpha$. These
calculations, using Laplace’s equation, pin the golden aspect ratio to be
around $\alpha^{\star}=1.07$ for the cylindrical box (and around 1.2 for a
rectangular box). Furthermore, these values of the golden aspect ratio are
remarkably close to the estimate in Eq. 16. We also want to examine more
realistic simulations, and thus we now examine the PNP equations, which are
widely employed to compute ion channel behavior Eisenberg (1996); Singer and
Norbury (2009); Zheng and Wei (2011); Burger (2011); Modi  _et al._ (2012);
Lin  _et al._ (2013) and MD simulations.

Figure 7:  The resistivity along the $z$-axis for various applied potentials
and cell dimensions $L$. The pore radius is $a=1$ nm and the KCl concentration
is $0.1$ mol/L ($\gamma_{b}=710$ M$\Omega\cdot$nm). Large applied voltages
enhance the density of charge carriers near the pore, decreasing the
resistivity. Since the access resistance contribution is largest near the
pore, see Eq. II, this enhancement lowers the overall access resistance and
gives an apparent depression to the golden aspect ratio when examining small
simulation cells. Figure 8:  PNP solution of the resistance versus $L$ at
applied potential of 1 V and KCl concentrations 0.1 mol/L (top panel) and 1.0
mol/L (bottom panel) for cylindrical cells with $h_{p}=1$ nm, $a=1$ nm, and
aspect ratios ($\alpha=$1.4, 1.1 and 0.8). The solid lines show the fit with
$f$ and $R_{\infty}$ as fitting parameters (shown in bold) and the dashed
lines show the extrapolation. Since the applied field is large, the resistance
at smaller $L$ ($<$ 100 nm) is highly nonlinear and cannot be scaled to
extract $R_{\infty}$ with Eq. 12 as is done for the smaller applied voltages
in Fig. 6. Hence, a larger $L$ is necessary to employ the scaling form when
the electric field is large and ion concentrations are small. Note that the
resistance for the 0.1 mol/L solution is substantially smaller than at lower
voltages. For 0.1 mol/L solution, the error of the fit for $f$ and
$R_{\infty}$ are about 3 % and 0.03 % respectively, whereas for 1.0 mol/L
solution the respective errors are less than 1 % and 0.01 %.

PNP Solution
We test Eqs. 12 for PNP simulations with two different concentrations, 0.1
mol/L and 1.0 mol/L of KCl solution, and applied voltages of 10 mV, 100 mV and
300 mV. Fig. 6 shows that for 10 mV and 100 mV we obtain a similar golden
aspect ratio as for the homogeneous medium solution. However, this special
ratio seems to decrease for the larger voltage of 300 mV. Clearly, there are
nonlinear effects coming into play when the voltage increases and/or the
concentration is insufficient to screen the field without substantially
perturbing the medium. We will show that this decrease is just an apparent
decrease due to non-scaling [in the context of Eq. 12] finite-size effects.
Larger $L$ will restore the scaling form.

We first want to understand the origin of these effects, which requires that
we consider the formation of concentration gradients Levadny  _et al._ (1998).
As seen in Fig. 7, the resistivity ($\gamma=1/\sum c_{\nu}\mu_{\nu}$) is
fairly constant along the axis of the simulation cell for small voltages and
thus we get the same result as in the homogeneous case. For higher voltages,
however, $\gamma$ is smaller near the pore as a result of classical Wien
effect Robinson and Stokes (2002); Onsager and Kim (1957). However, the
scaling form, Eq. 12, can still be employed as long as the resistivities of
different regions do not change drastically. Relaxing the homogeneity
assumption that went into the derivation of Eq. 12, we set $\gamma_{b}$,
$\gamma_{a}$, and $\gamma_{p}$ to be the resistivity in the bulk-like region,
the access-like region, and in the pore, respectively. Then, the total
resistance (for $H\gtrsim L$) is

$\displaystyle R=$ $\displaystyle
R_{\mathrm{access}}+R_{\mathrm{bulk}}+R_{\mathrm{intermediate}}+R_{\mathrm{pore}}$
$\displaystyle=$ $\displaystyle
2\left(\frac{\gamma_{a}}{4a}-\frac{\gamma_{a}}{\pi
f_{1}L}\right)+\gamma_{b}\left(\frac{\alpha
L-f_{2}L}{\mathcal{G}L^{2}}\right)+\frac{\gamma_{b}f_{3}}{\mathcal{G}L}+\frac{\gamma_{p}h_{p}}{\pi
a^{2}}$ $\displaystyle=$
$\displaystyle\frac{\gamma_{b}}{\mathcal{G}L}\left[\alpha-f\left(\frac{\gamma_{a}}{\gamma_{b}}\right)\right]+R^{\prime}_{\infty},$
(19)

where $R^{\prime}_{\infty}=\gamma_{a}/2a+\gamma_{p}h_{p}/\pi a^{2}$ and

$f\left(\frac{\gamma_{a}}{\gamma_{b}}\right)=\frac{2\mathcal{G}\gamma_{a}}{\pi
f_{1}\gamma_{b}}+f_{2}-f_{3}.$ (20)

Within this simplified – “compartmentalized” – inhomogeneous system, the
resistance of the infinite system and the golden aspect ratio are linear
function of $\gamma_{a}$, decreasing as $\gamma_{a}$ decreases. This
qualitatively explains why we see the apparent decrease in the golden aspect
ratio – it is the depression of the access contribution that decreases the
required height of the bulk-like region that “balances” it. Of course, the
resistivity is not constant throughout the access-like region and, at very
high voltage and low concentration, is not even in the bulk-like region. Thus,
this simplified inhomogeneous model will not accurately determine the apparent
golden aspect ratio. However, it does not need to: The extensive decrease in
resistivity is due to the fact that the simulation cells that are really small
can not sustain large potential drops without a global perturbation to the
medium. Examining much larger cells shows that the scaling form is restored,
see Fig. 8, which is due to the fact that the perturbation to the medium
becomes localized (in dimensionless terms) around the pore, see the curves for
different $L$ in Fig. 7.

In other words, there are finite-size effects and there are finite-size
effects: The perturbation to the medium is a non-scaling finite size effect
and the simulation cells need to be sufficiently large to accommodate the
perturbed region. An alternative scaling ansatz could potentially be developed
to handle such cases, e.g., including the perturbing effect of local fields.
Its applicability would be more limited (unlike the general forms, Eq. 12 or
Eq. 8, which, however, may require larger $L$ to reach the scaling regime) but
still helpful. Finite-size effects related to truncating the access region or
asymmetric dimensions can be removed via Eq. 12 or Eq. 8 (for $H\gtrsim L/f$),
respectively.

MD Solution
We test Eqs. 12 for MD simulations with 1 mol/L KCl, a 1.81 nm radius graphene
pore, and a 1 V applied bias. Figure 9 shows results for three different
aspect ratios, 1.9,1.2 and 0.5, and a cross-sectional length ranging from 9.6
nm to 14.4 nm. Both equations hold for these all-atom simulations. Moreover,
we obtain the same golden aspect ratio for MD simulation as we do for the
continuum calculations with a rectangular simulation cell.

Figure 9:  MD results of the resistance versus $L$ at an applied bias of 1 V
and KCl concentration 1.0 mol/L, $a=1.81$ nm, and aspect ratios 1.9, 1.2 and
0.5 (from top to bottom). The raw resistance from MD is shown as open circles.
Due to the equilibration, the final cell aspect ratios are not exactly 1.9,
1.2, and 0.5. The slight change in height can be corrected for in order to set
$H=\alpha L+h_{p}$. The data with error bars show this slightly corrected
data. The error bars represent $\pm$ 1 block standard error, see Ref.
Grossfield and Zuckerman (2009), which reflects error due to finite sampling
time. The legends show the fits to Eq. (12) with the fit parameters shown in
bold. The golden aspect ratio for MD is $\approx 1.2$, which is the same as
that from continuum simulations (for a rectangular cell). The scaling and/or
golden aspect ratio simulations yield a resistance in agreement with the
classical expression, $2R_{\mathrm{Hall}}+R_{\mathrm{pore}}$, as we
demonstrated with the scaling procedure only in Sahu et. al. Sahu and Zwolak
(2018).

##  IV DISCUSSION

Ion transport properties can be difficult to compute with all-atom molecular
dynamics. In the dehydration limit, for instance, large free energy barriers
exist, which entails very long (microseconds or more) simulations to
accumulate enough ion crossing events to determine the current with reasonable
accuracy Sahu  _et al._ (2017); Sahu and Zwolak (2017). This applies even at
elevated voltages [simulations are often done at 3 V or higher Sathe  _et al._
(2011)], which risks introducing nonlinearities and prohibits a direct
comparison with experiment. Contextual properties of ion transport –
fluctuations of the pore or membrane, localized dipoles or charges, van der
Waals interactions and dehydration, geometrically “imperfect” pore shapes,
etc. – make these simulations even more difficult to understand and control
(e.g., slow fluctuations of pore structure require long times to acquire the
requisite statistics). To capture these effects, therefore, one would like to
keep the simulation size as small as possible but without introducing
excessive errors into the computation.

We have developed a scaling theory to extract both the pore and access
contributions to resistance without going to excessively large simulation
cells (i.e., $L<16$ nm for linear response). This theory was originally
employed to determine the access contribution – including how to properly
define the pore size in the presence of contextual properties – to graphene
nanopores with all-atom MD Sahu and Zwolak (2018). While this does not correct
for force-field effects and related issues (sampling, conformational changes,
etc.), it does solve for the geometry/dimensions of the setup. Using this
theory, we demonstrated that the “golden aspect ratio” – a special aspect
ratio where finite size effects are eliminated – exists in both continuum and
all-atom simulations (and we expect in Brownian dynamics as well). Previously,
we only used an aspect ratio close to an initial estimate of the golden aspect
ratio (and, indeed, we had little variation of $R$ with $L$) but large enough
to ensure we decrease from above to $R_{\infty}$ Sahu and Zwolak (2018).

The golden aspect ratio gives a completely flat (i.e., without non-monotonic
or oscillatory behavior) resistance versus simulation cell size. The value for
the golden aspect ratio is about $\alpha^{\star}=1.07$ for cylindrical cells
and 1.2 for rectangular (for both continuum and all-atom simulations). It is
relatively insensitive to voltages and ion concentrations within
experimentally relevant regimes (i.e., around 0.1 V applied voltages and 0.1
mol/L to 1.0 mol/L concentrations). For large voltages and small
concentrations, simulation cells must be sufficiently large to remove non-
scaling finite-size effects, after which the scaling and the golden aspect
ratio are restored.

This scaling approach will be most useful when examining pores that do not
have very high pore resistances – due to dehydration or a long length channel
– so that the access resistance becomes a dominant or equal contributor rather
than a correction to the resistance. Graphene and other atomically thin pores
have a dominant or substantial access contribution all the way down to the
dehydration limit Sahu and Zwolak (2018), and thus the scaling approach is
essential. Biological channels can fall into either category, sometimes
requiring a correct determination of the access contribution. For instance,
$\alpha$-hemolysin pores have a 1 G$\Omega$ resistance (at 1 mol/L KCl)
Aksimentiev and Schulten (2005), compared to an estimated [from Eq. 1] 40
M$\Omega$ access contribution. On the other hand, a sodium channel of radius
0.3 nm and effective length 0.5 nm has both an access and pore resistance that
are nearly equal Hille (1968). In low salt concentrations, access resistance
may become the dominating resistance in a wide variety of cases Alcaraz  _et
al._ (2017). While we have not simulated a biological ion channel here, a
direct application of the golden aspect ratio would take the simulation cell
height to be $1.2\cdot L+h_{p}$, with $h_{p}$ the channel length. For
instance, for $\alpha$-hemolysin as a model biological channel, $h_{p}\approx
10$ nm, the simulation cell height should be roughly on the order of 27 nm
(for the cross-sectional dimension to be reasonably larger than channel width
and protein assembly, $L$ should be about 14 nm Aksimentiev and Schulten
(2005)). However, this idea requires testing with different biological
channels, as in many cases (including $\alpha$-hemolysin), the membrane
thickness and channel length are not the same. This distorts the equipotential
surfaces. Given the ease with which the scaling approach and golden aspect
ratio can be employed (e.g., setting $H\approx 1.2\cdot L+h_{p}$ and
simulating a handful of different $L$), it is reasonable to use it first even
in the cases the where access resistance is a correction.

Care must be taken, though, to work with sufficiently large $L$ in
simulations, especially all-atom MD. For instance, a simulation cell that is
too small may not have enough ions to properly screen localized dipoles or
charges (or otherwise give a cell boundary too close to these local
electrostatic variations, preventing proper Debye screening). Similarly, the
partitioning of the voltage drop across the cell can yield a strong field
across the pore for small simulation cells (as there is less voltage drop in
the bulk medium), which can introduce non-linearities in the transport
properties. This further supports the use of the scaling analysis, as it will
also allow one to identify these “non-scaling” finite size effects and
eliminate them when going to large enough simulation cells. Moreover, it goes
without saying that computations should be done as close to experimental and
biologically-relevant conditions as possible. Our work gives new reasons to do
simulations in the linear response limit, as it is only then that some non-
scaling finite size effects are negligible for simulation cell sizes in the 16
nm regime.

In addition, the classical form of access resistance assumes an infinite bulk,
which fails for micro- and nano-scale systems. In particular, if one has a
narrow (micro- or nano-fluidic) constriction leading up to a membrane/pore or
some other “active” region, then access resistance no longer follows Hall’s
expression, but rather has to be corrected for the geometric setup Green  _et
al._ (2016). Our approach gives a method to map simulations to experimental
geometries (or vice versa), thus allowing for contextual properties to be
simulated with reasonable computational power.

The primary power of MD is to study contextual properties of pores, such as
the effect of atomic/molecular fluctuations of the pore/membrane and the
presence of (partial) charges, but requires extensive computational resources.
The golden aspect ratio and scaling approaches will allow the quantitative
extraction of both pore and access contributions to the resistance with
minimal resources. Its use has already shed light on how to define pore size
for geometrically imperfect pores, including edge fluctuations and
dehydration/van der Waals interactions Sahu and Zwolak (2018). By extension,
this will give the first calculated values of access resistance in the
presence of, e.g., molecular-scale/protein fluctuations in biological ion
channels and other contextual properties (charges/dipoles) and open a new era
in comparing computed and measured values of resistance.

*Appendix: Access resistance in a finite region

Here, we calculate access resistance for a finite size region while keeping
the ellipsoidal symmetry. This problem can be solved using the rotational
elliptic coordinates. We begin with oblate spheroidal coordinates, ($\mu$,
$\nu$, $\phi$), which are defined in terms of the Cartesian coordiates as

$\displaystyle x$ $\displaystyle=a\cosh\mu\,\cos\nu\,\cos\phi$ $\displaystyle
y$ $\displaystyle=a\cosh\mu\,\cos\nu\,\sin\phi$ $\displaystyle z$
$\displaystyle=a\sinh\mu\,\sin\nu\,$ (21)

Rotational elliptic coordinates, $\xi=\sinh\mu$ and $\eta=\sin\nu$, are
related to the cylindrical coordinates as

$\displaystyle z=a\xi\eta$
$\displaystyle\rho=a\sqrt{(1+\xi^{2})(1-\eta^{2})},$ (22)

where $a$ is the radius of the pore. Laplace’s equation in this coordinate
system is

$\frac{\partial}{\partial\xi}\left[(1+\xi^{2})\frac{\partial
V}{\partial\xi}\right]+\frac{\partial}{\partial\eta}\left[(1-\eta^{2})\frac{\partial
V}{\partial\eta}\right]=0$ (23)

This equation can be solved via separation of variables, $V=H(\eta)\Xi(\xi)$,
which gives

$\frac{\partial}{\partial\xi}\left[(1+\xi^{2})\frac{\partial\Xi}{\partial\xi}\right]+\lambda\,\Xi=0$
(24)

and

$\frac{\partial}{\partial\eta}\left[(1-\eta^{2})\frac{\partial
H}{\partial\eta}\right]-\lambda\,H=0,$ (25)

where $\lambda$ is a constant. The boundary condition (ii), $\partial
V/\partial\eta=0$ at $\eta=0$, yields $\lambda=0$. Integrating Eq. 24 thus
gives

$V=H\,\Xi=H\,\int\frac{b_{1}}{1+\xi^{2}}d\xi=H\,b_{1}\tan^{-1}\xi+b_{2},$ (26)

where $b_{1}$ and $b_{2}$ are the constants of integration. Using the boundary
condition (i), $V(\xi=0)=0$, we get $b_{2}=0$. This simplifies the potential
to $V=H\,b_{1}\tan^{-1}\xi$. The boundary condition (iii), $V(\xi=l/a)=V_{0}$,
gives $V_{0}=H\,b_{1}\tan^{-1}(l/a)$. Taking the ratio of these two
expressions gives Eq. 4 of the main text,

$\frac{V}{V_{0}}=\frac{\tan^{-1}\xi}{\tan^{-1}(l/a)}.$ (27)

The current through the pore is

$\displaystyle I=$ $\displaystyle\frac{2\pi}{\gamma}\int_{0}^{a}\frac{\partial
V}{\partial z}\Big{|}_{z=0}\rho
d\rho=\frac{2\pi}{\gamma}\int_{0}^{a}\frac{1}{a\eta}\frac{\partial
V}{\partial\xi}\Big{|}_{\xi=0}\ \rho d\rho$ $\displaystyle=$
$\displaystyle\frac{2\pi aV_{o}}{\gamma\tan^{-1}(l/a)}.$ (28)

Equations 27 and 28 give Eq. II.

$\displaystyle R_{\mathrm{access/2}}$
$\displaystyle=\frac{\gamma\tan^{-1}(l/a)}{2\pi a}\approx
R_{\mathrm{Hall}}-\frac{\gamma}{2\pi l}.$

## Acknowledgments

We thank J. Elenewski, D. Gruss, and C. Rohmann for helpful comments. S. S.
acknowledges support under the Cooperative Research Agreement between the
University of Maryland and the National Institute of Standards and Technology
Center for Nanoscale Science and Technology, Award 70NANB14H209, through the
University of Maryland.

## References

  * Hille (2001) B. Hille,  _Ion channels of excitable membranes_ , Vol. 507 (Sinauer Sunderland, MA, 2001).
  * Kandel  _et al._ (2000) E. R. Kandel, J. H. Schwartz, T. M. Jessell, S. A. Siegelbaum, A. J. Hudspeth, _et al._ ,  _Principles of neural science_ , Vol. 4 (McGraw-hill New York, 2000).
  * Doyle  _et al._ (1998) D. A. Doyle, J. a. M. Cabral, R. A. Pfuetzner, A. Kuo, J. M. Gulbis, S. L. Cohen, B. T. Chait, and R. MacKinnon, Science 280, 69 (1998).
  * Zwolak and Di Ventra (2008) M. Zwolak and M. Di Ventra, Rev. Mod. Phys. 80, 141 (2008).
  * Heerema and Dekker (2016) S. J. Heerema and C. Dekker, Nat. Nanotechnol. 11, 127 (2016).
  * Baker (2000) R. W. Baker,  _Membrane Technology_ (Wiley Online Library, 2000).
  * Mulder (2012) M. Mulder,  _Basic principles of membrane technology_ (Springer Science & Business Media, 2012).
  * Garaj  _et al._ (2010) S. Garaj, W. Hubbard, A. Reina, J. Kong, D. Branton, and J. Golovchenko, Nature 467, 190 (2010).
  * Merchant  _et al._ (2010) C. A. Merchant, K. Healy, M. Wanunu, V. Ray, N. Peterman, J. Bartel, M. D. Fischbein, K. Venta, Z. Luo, A. T. C. Johnson, and M. Drndić, Nano Lett. 10, 2915 (2010).
  * Schneider  _et al._ (2010) G. F. Schneider, S. W. Kowalczyk, V. E. Calado, G. Pandraud, H. W. Zandbergen, L. M. Vandersypen, and C. Dekker, Nano Lett. 10, 3163 (2010).
  * Heiranian  _et al._ (2015) M. Heiranian, A. B. Farimani, and N. R. Aluru, Nat. Commun. 6, 8616 (2015).
  * Walker  _et al._ (2017) M. I. Walker, K. Ubych, V. Saraswat, E. A. Chalklen, P. Braeuninger-Weimer, S. Caneva, R. S. Weatherup, S. Hofmann, and U. F. Keyser, ACS Nano 11, 1340 (2017).
  * Joshi  _et al._ (2014) R. Joshi, P. Carbone, F.-C. Wang, V. G. Kravets, Y. Su, I. V. Grigorieva, H. Wu, A. K. Geim, and R. R. Nair, Science 343, 752 (2014).
  * Abraham  _et al._ (2017) J. Abraham, K. S. Vasu, C. D. Williams, K. Gopinadhan, Y. Su, C. T. Cherian, J. Dix, E. Prestat, S. J. Haigh, I. V. Grigorieva, A. K. Geim, and R. R. Nair, Nat. Nanotechnol. 12, 546 (2017).
  * Esfandiar  _et al._ (2017) A. Esfandiar, B. Radha, F. C. Wang, Q. Yang, S. Hu, S. Garaj, R. R. Nair, A. K. Geim, and K. Gopinadhan, Science 358, 511 (2017).
  * Cohen-Tanugi and Grossman (2012) D. Cohen-Tanugi and J. C. Grossman, Nano Lett. 12, 3602 (2012).
  * Surwade  _et al._ (2015) S. P. Surwade, S. N. Smirnov, I. V. Vlassiouk, R. R. Unocic, G. M. Veith, S. Dai, and S. M. Mahurin, Nat. Nanotechnol. 10, 459 (2015).
  * Jiang  _et al._ (2009) D. Jiang, V. R. Cooper, and S. Dai, Nano Lett. 9, 4019 (2009).
  * Kim  _et al._ (2013) H. W. Kim, H. W. Yoon, S.-M. Yoon, B. M. Yoo, B. K. Ahn, Y. H. Cho, H. J. Shin, H. Yang, U. Paik, S. Kwon, J.-Y. Choi, and H. B. Park, Science 342, 91 (2013).
  * Clarke  _et al._ (2009) J. Clarke, H.-C. Wu, L. Jayasinghe, A. Patel, S. Reid, and H. Bayley, Nat. Nanotechnol. 4, 265 (2009).
  * Si and Aksimentiev (2017) W. Si and A. Aksimentiev, ACS Nano 11, 7091 (2017).
  * Zwolak  _et al._ (2009) M. Zwolak, J. Lagerqvist, and M. Di Ventra, Phys. Rev. Lett. 103, 128102 (2009).
  * Zwolak  _et al._ (2010) M. Zwolak, J. Wilson, and M. Di Ventra, J. Phys.: Condens. Matter 22, 454126 (2010).
  * Sahu  _et al._ (2017) S. Sahu, M. Di Ventra, and M. Zwolak, Nano Lett. 17, 4719 (2017).
  * Sahu and Zwolak (2017) S. Sahu and M. Zwolak, Nanoscale 9, 11424 (2017).
  * Hall (1975) J. E. Hall, J. Gen. Physiol. 66, 531 (1975).
  * Hille (1968) B. Hille, J. Gen. Phsiol. 51, 199 (1968).
  * Läuger (1976) P. Läuger, Biochim. Biophys. Acta 455, 493 (1976).
  * Alcaraz  _et al._ (2017) A. Alcaraz, M. L. López, M. Queralt-Martín, and V. M. Aguilella, ACS Nano 11, 10392 (2017).
  * Vodyanoy and Bezrukov (1992) I. Vodyanoy and S. M. Bezrukov, Biophys. J. 62, 10 (1992).
  * Bezrukov and Vodyanoy (1993) S. M. Bezrukov and I. Vodyanoy, Biophys. J. 64, 16 (1993).
  * Bezrukov  _et al._ (1996) S. M. Bezrukov, I. Vodyanoy, R. A. Brutyan, and J. J. Kasianowicz, Macromolecules 29, 8517 (1996).
  * Maffeo  _et al._ (2012) C. Maffeo, S. Bhattacharya, J. Yoo, D. Wells, and A. Aksimentiev, Chem. Rev. 112, 6250 (2012).
  * Sahu and Zwolak (2018) S. Sahu and M. Zwolak, Phys. Chem. Chem. Phys. 20, 4646 (2018).
  * Eisenberg (2010) B. Eisenberg, J. Phys. Chem. C 114, 20719 (2010).
  * Yoo and Aksimentiev (2015) J. Yoo and A. Aksimentiev, J Phys Chem Lett 6, 4680 (2015).
  * Maxwell (1881) J. C. Maxwell,  _A treatise on electricity and magnetism_ , Vol. 1 (Clarendon press, 1881).
  * Gray and Mathews (1895) A. Gray and G. B. Mathews,  _A treatise on Bessel functions and their applications to physics_ (Macmillan and Company, 1895).
  * Brown and Escombe (1900) H. T. Brown and F. Escombe, Proc. Roy. Soc. London 67, 124 (1900).
  * Gröber (1921) H. Gröber,  _Die Grundgesetze der Wärmeleitung und des Wärmeüberganges: ein Lehrbuch für Praxis und technische Forschung_ (Springer-Verlag, Berlin, 1921).
  * Holm (1958) R. Holm,  _The contact resistance. General theory_ (Springer, 1958).
  * Newman (1966) J. Newman, J. Electrochem. Soc. 113, 501 (1966).
  * Sharvin (1965) Y. V. Sharvin, J. Exp. Theor. Phys. 21, 655 (1965).
  * Wexler (1966) G. Wexler, Proc. Phys. Soc. 89, 927 (1966).
  * Nikolić and Allen (1999) B. Nikolić and P. B. Allen, Phys. Rev. B 60, 3963 (1999).
  * Aguilella-Arzo  _et al._ (2005) M. Aguilella-Arzo, V. M. Aguilella, and R. S. Eisenberg, Eur. Biophys. J. 34, 314 (2005).
  * Luchinsky  _et al._ (2009) D. Luchinsky, R. Tindjong, I. Kaufman, P. McClintock, and R. Eisenberg, Phys. Rev. E 80, 021925 (2009).
  * Peskoff and Bers (1988) A. Peskoff and D. Bers, Biophys. J. 53, 863 (1988).
  * Panday  _et al._ (2016) N. Panday, G. Qian, X. Wang, S. Chang, P. Pandey, and J. He, ACS Nano 10, 11237 (2016).
  * Read  _et al._ (2009) M. B. Read, J. H. Lang, A. H. Slocum, and R. Martens, in  _Proceedings of the 55th IEEE Holm Conference on Electrical Contacts_ (IEEE, 2009) pp. 303–309.
  * Braunovic  _et al._ (2006) M. Braunovic, N. K. Myshkin, and V. V. Konchits,  _Electrical contacts: Fundamentals, applications and technology_ (CRC press, 2006).
  * Abramowitz and Stegun (1964) M. Abramowitz and I. A. Stegun, eds.,  _Handbook of mathematical functions: with formulas, graphs, and mathematical tables_ (Courier Corporation, 1964) pp. 752–770.
  * Mikic (1967) B. B. Mikic,  _Thermal contact resistance._ , Ph.D. thesis, Massachusetts Institute of Technology (1967).
  * Cooper  _et al._ (1969) M. Cooper, B. Mikic, and M. Yovanovich, Int J Heat Mass Transfer 12, 279 (1969).
  * Bhattacharya  _et al._ (2011) S. Bhattacharya, J. Muzard, L. Payet, J. Mathé, U. Bockelmann, A. Aksimentiev, and V. Viasnoff, J. Phys. Chem. C 115, 4255 (2011).
  * Phillips  _et al._ (2005) J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kale, and K. Schulten, J. Comput. Chem. 26, 1781 (2005).
  * Feller and MacKerell (2000) S. E. Feller and A. D. MacKerell, J. Phys. Chem. B 104, 7510 (2000).
  * Jorgensen  _et al._ (1983) W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983).
  * Eisenberg (1996) R. Eisenberg, J. Membr. Biol. 150, 1 (1996).
  * Singer and Norbury (2009) A. Singer and J. Norbury, SIAM J. Appl. Math. 70, 949 (2009).
  * Zheng and Wei (2011) Q. Zheng and G.-W. Wei, J. Chem. Phys. 134, 194101 (2011).
  * Burger (2011) M. Burger, Inverse Prob. 27, 083001 (2011).
  * Modi  _et al._ (2012) N. Modi, M. Winterhalter, and U. Kleinekathoefer, Nanoscale 4, 6166 (2012).
  * Lin  _et al._ (2013) G. Lin, W. Liu, Y. Yi, and M. Zhang, SIAM J. Appl. Dyn. Syst. 12, 1613 (2013).
  * Levadny  _et al._ (1998) V. Levadny, V. M. Aguilella, and M. Belaya, Biochimica et Biophysica Acta (BBA) - Biomembranes 1368, 338 (1998).
  * Robinson and Stokes (2002) R. A. Robinson and R. H. Stokes,  _Electrolyte solutions_ (Courier Corporation, 2002).
  * Onsager and Kim (1957) L. Onsager and S. K. Kim, J. Phys. Chem. 61, 198 (1957).
  * Grossfield and Zuckerman (2009) A. Grossfield and D. M. Zuckerman, Annu. Rep. Comput. Chem. 5, 23 (2009).
  * Sathe  _et al._ (2011) C. Sathe, X. Zou, J.-P. Leburton, and K. Schulten, ACS Nano 5, 8842 (2011).
  * Aksimentiev and Schulten (2005) A. Aksimentiev and K. Schulten, Biophys. J. 88, 3745 (2005).
  * Green  _et al._ (2016) Y. Green, R. Eshel, S. Park, and G. Yossifon, Nano Lett. 16, 2744 (2016).