Table Of ContentMultiscale modeling of a rectifying bipolar nanopore: Comparing
Poisson-Nernst-Planck to Monte Carlo
Bartl omiejMatejczyk1, Mo´nika Valisk´o2, Marie-ThereseWolfram1,3, Jan-FrederikPietschmann4 andDezs˝o Boda2,5∗
1Department of Mathematics, University of Warwick, CV4 7AL Coventry, United Kingdom
2Department of Physical Chemistry, University of Pannonia, P.O. Box 158, H-8201 Veszpr´em, Hungary
3Radon Institute for Computational and Applied Mathematics,
Austrian Academy of Sciences, Altenberger Strasse 69, 4040 Linz, Austria
4Institute for Computational and Applied Mathematics, WWU Mu¨nster, Mu¨nster 48149, Germany and
5Institute of Advanced Studies K˝oszeg (iASK), Chernel st. 14. H-9730 K˝oszeg, Hungary
7 (Dated: January 24, 2017)
1
In the framework of a multiscale modeling approach, we present a systematic study of a bipolar
0
rectifying nanopore using a continuum and a particle simulation method. The common ground in
2
thetwomethodsistheapplication oftheNernst-Planck(NP)equationtocomputeiontransportin
n theframeworkoftheimplicit-waterelectrolytemodel. ThedifferenceisthatthePoisson-Boltzmann
a theory is used in the Poisson-Nernst-Planck (PNP) approach, while the Local Equilibrium Monte
J Carlo (LEMC) method is used in the particle simulation approach (NP+LEMC) to relate the con-
3 centration profile to the electrochemical potential profile. Since we consider a bipolar pore which
2 is short and narrow, we perform simulations using two-dimensional PNP. In addition, results of a
non-linear version of PNP that takes crowding of ions into account are shown. We observe that
] the mean field approximation applied in PNP is appropriate to reproduce the basic behavior of
t
f the bipolar nanopore (e.g., rectification) for varying parameters of the system (voltage, surface
o
charge, electrolyte concentration, and pore radius). We present current data that characterize the
s
nanopore’s behavior as a device, as well as concentration, electrical potential, and electrochemical
.
t potential profiles.
a
m
PACSnumbers: 87.16.dp,02.70.-c,05.10.Ln,07.05.Tp
-
d
n
I. INTRODUCTION µEX(r) term is the excess chemical potential that de-
o i
scribes all the interactions acting between the particles
c
[ In this paper, we compare Poisson-Nernst-Planck formingthesystemandalltheinteractionswithexternal
(PNP) theory with particle simulations for ionic trans- forces (including an applied electrical potential). PNP
1
v port through a rectifying bipolar nanopore. Both meth- defines the excess term as the interaction with the mean
4 ods use the Nernst-Planck (NP) transport equation to electricfieldproducedbyallthefreechargesandinduced
4 describe the ionic flux of i={1,2} species: charges. Thus the electrochemical potential in the case
3 of PNP is
6 1
0 ji(r)=−kTDi(r)ci(r)∇µi(r), (1) µPiNP(r)=µ0i +kT lnci(r)+zieΦ(r), (3)
.
1 whereji(r)istheparticlefluxdensityofionicspeciesi,k where zi is the ionic valence, e the elementary charge,
0 theBoltzmann’sconstant,T thetemperature,andD (r) andΦ(r)thetotalmeanpotential. Themissingtermcan
i
7 the diffusion coefficient profile. The main difference be- be identified with what is beyond meanfield (BMF) and
1
tween the two techniques is that PNP makes use of the quantifiesthedifferencebetweenPNPandasolutionthat
:
v Poisson-Boltzmann(PB) theory to relate the concentra- isaccuratefromthepointofviewofstatisticalmechanics:
Xi tion profile, ci(r), to the electrochemical potential pro- µ (r)=µPNP(r)+µBMF(r). (4)
file, µ (r), while the particle simulation method uses the i i i
r i
a LocalEquilibriumMonte Carlo(LEMC)technique1–4 to In the implicit solvent framework used here the BMF
establish this relation. The particle simulation method term includes the volume exclusion effects (hard sphere
includes all the ionic correlations that are beyond the effects) due to the finite size ofthe ions andelectrostatic
mean field approximation applied in PNP. The differ- correlations that are beyond the mean-field level. This
ence between the two approaches can be quantified by partitioning has been used to study selective adsorption
considering the electrochemical potential of ions at electrodes5 and in ion channels6–8.
It is also usual to break the electrochemical poten-
µ (r)=µ0+kT lnc (r)+µEX(r), (2)
i i i i tialintoa chemicalandanelectricalcomponentthatare
where µ0 is a standard chemical potential, a constant loosely identified with the chemical and electrical works
i
term that does not appear in the calculations. The needed to bring an ion from one medium to the other:
µ (r)=µCH(r)+µEL(r), (5)
i i i
where the EL term can be identified with z eΦ(r), while
i
∗Authorforcorrespondence: [email protected] the CH term can be identified with µ0i + kT lnci(r) +
2
µBMF(r). Although these two terms cannot be sepa- could be termed as NP+PB, but we stay with the usual
i
rated in experiments9–11, the separation is possible in name,PNP.Poisson’sequationis satisfiedinbothmeth-
computational studies because Φ(r) can be determined. ods. InPNP,itissolvedineveryiteration,whileitisau-
In PNP, where µBMF(r) = 0, the CH term is just tomatically fulfilled in LEMC because Coulomb’s law is
i
µID(r)=µ0+kTlnc (r), the ideal expression (ID). usedtohandleelectrostaticsinthesimulations(including
i i i
The PNP model is based on an approximate mean theappliedfieldinthe frameworkoftheInducedCharge
fieldapproachwithalltheadvantagesanddisadvantages. Computationmethod27,28). Bothapproachesprovideap-
First of all, the mean field method does not consider the proximate indirect solutions for the dynamical problem
particlesasindividualentities, butworkswiththeir con- through the NP equation. Direct simulation of ionic
centrationprofiles whichcan be understoodas the prob- transportin the implicit solventframeworkis commonly
ability of finding an ion at a specific point in space and donebyBrownianDynamics(BD)simulations29–31. The
time. Thisprobabilitydependsontheinteractionenergy main difference between NP+LEMC and PNP is the
of the ion with the average (mean) electrical potential way they handle the statistical mechanical problem of
produced by all the charges in the system, including all establishing the closure between c (r) and µ (r). The
i i
theions. Two-andmany-bodycorrelationsbetweenions, NP+LEMC technique provides a solution on the basis
therefore, are neglected in PNP. The ions are treated as of particle simulations that contain all the correlations
point charges omitting their size. ignored by PNP. The main goal of our study is to dis-
In this work, we also use a non-linear variant of PNP cuss the effects of the approximations applied in PNP
(denotedbynPNPfromnowon)thatcanbederived(for- for different sets of physical parameters. Comparing to
mally) from a discrete hopping model12. Similar models NP+LEMC results makes it possible to focus on the ap-
have been derived by Bikerman13 and Li14. In all these proximationsappliedinthestatisticalmechanicalpartof
models Eqs. (1) and (3) are replaced by thePNPtheory(thePBtheory),becauseNPiscommon
in them. If we want to reveal the magnitude and nature
1
jnPNP(r)=− D (r)c (r)α(r)∇µnPNP(r) (6) of errors resulting from the application of the approxi-
i kT i i i mative NP equation instead of simulating ion transport
directly,weneedtocomparetoBD32,33. Comparisonbe-
and
tween BD and NP+LEMC results will be published in a
µnPNP(r)=µPNP(r)−kT lnα(r), (7) separate paper.
i i
We apply our methods to a bipolar nanopore that is a
where suitable case study for our purpose. Bipolar nanopores
have an asymmetrical surface charge distribution on the
c (r) c (r)
α(r)=1− 1 − 2 , (8) pore wall changing sign along the central axis of the
c c
max max pore. Pore regions with opposite surface charges can be
achieved by chemical modifications. For example, in the
andwehavechosenc =61.5mol/dm3 asamaximum
max
case of PET nanopores, carboxyl groups can be trans-
value for the concentation at close packing. The scaling
formed into amino groups by a coupling agent34. The
factor, α(r) approaches 1 as c →0, so nPNP turns into
i
surface potential can also be regulated similarly to field-
PNP in this limit. Another approach to overcome these
effecttransistorsiftheporewallsaremadeofconducting
limitation is density functional theory (DFT)15–23 which
materials.
includesadditionaltermsinthefreeenergythattakecare
The reason of choosing the bipolar nanopore for this
of the interactions in the BMF term. A recent review
discusses different possibilities to account for ions size24. comparative work is that the source of rectification in
this case is purely electrostatic in nature and thus a ro-
In most of the literature, a one-dimensional (1D) re-
bust effect. Therefore, we can afford a short nanopore
ductionofPNPisused. Thisgivesgoodresultsespecially
forlongandnarrownanopores25asitsderivationisbased (only 6 nm in length) that can be handled with LEMC.
Inthecaseofconicalnanopores,whereonlyageometrical
ontheassumptionthattheradiusissignificantlysmaller
asymmetry is present, long pores are needed to produce
than the length. Its advantage is that it requires less
aconsiderableeffectwhichmakesitcomputationallyun-
computationaleffort and makes the computation of long
feasible.
pores possible. In this work, however, we use a two di-
mensional(2D)PNP(respectivelynPNP)modelwhichis Although bipolar nanopores have been studied exten-
a suitable approximation to the three-dimensional (3D), sively using PNP25,34–46, we are not awareof any paper,
but rotationally symmetric, system studied here. Our where a direct comparison to particle simulations is dis-
model,furthermore,includesthebulkregionsandtheac- cussedforthissystem. Furthermore,mostofthoseworks
cessregionsattheentrancesofthenanopore,asopposed use1DPNP,whilewereportresultsof2DPNPhereand
tootherstudies25,26. SolvingPNPinalargerdomainand the comparison with nPNP is also completely new.
also in the radial dimension gives more accurate results. Particle simulations are necessary for narrow pores,
Summarised, we can couple the NP equation either to where ions are crowded and their size and the correla-
LEMC simulations or to the PB theory. The former is tions between them (the BMF term) matter. This is the
referredtoasthe NP+LEMCtechnique, while the latter case in ion channels, where the ions correlate strongly
3
with each other and with the charged amino acids along CoulombpotentialandascoefficientinthePoissonequa-
the ionic pathway. Although nanopores are larger in re- tion,respectively),whileitsdynamiceffectisincludedin
ality, the electrical double layers formed by the ions at the diffusioncoefficientinthe NPequation(D (r)inEq.
i
the pore walls overlap if the the Debye length is larger 1). The ions are point charges in PNP, while they are
than the pore radius. This occurs if the pore is narrow hard spheres (of radius 0.3 nm for both ions) with point
enough(suchasconicalnanoporesattheir tips) orif the charges at their centers in LEMC.
electrolyte is dilute. The nanopore is a cylinder of 6 nm in length with a
Thisworkbelongstoaseriesofstudies,whereweapply varying radius (R = 0.5−3 nm). It penetrates a mem-
a multiscale modeling approach47. We can create differ- brane that separates two bulk electrolytes. The walls of
ent models (with less or more details) and we can study the pore and the membrane are hard impenetrable sur-
these models with computational methods that fit the faces in the LEMC simulations (Fig. 1A), while they are
model. In another recent work,48, we compared results partoftheboundariesofthesolutiondomaininthePNP
of molecular dynamics (MD) simulations performed for calculations (Fig. 1B).
anall-atommodelincludingexplicitwatertoNP+LEMC Thediffusioncoefficientofthe ionsareusuallysmaller
calculations performed for the implicit-water model (the inside the pore than outside in the bulk regions. This
same model studied in this paper). We concluded that, finding was confirmed by our other study that com-
despite all the simplifications, the implicit-water model pares MD and NP+LEMC results48. Here, for sim-
provides an appropriate framework to study nanopores. plicity, we assigned Dbulk = 1.333 × 10−9 m2s−1 and
i
The link betweenthe twomodelinglevelsis the diffusion Dpore = 1.333×10−10 m2s−1 values in the bulk and in
i
coefficientprofile,Di(r),usedinNP+LEMCasaninput, the pore, respectively, for both ions.
whiletheMDsimulationscanprovideinformationabout
The chargeson the cylinder’s surface are partialpoint
this profile.
charges in the case of LEMC that are placed on grid
The advantage of NP+LEMC over MD is that it is
points whose average distance is about 0.25 nm. The
faster and can handle larger systems that are closer to
values of the partial charges depend on the prescribed
realistic length scales of nanodevices. From this point of
surface charge density, σ. The surface charge densities
view, PNP is even more advantageous, because it does
are included in PNP throughNeumann boundary condi-
not involve particle simulations, therefore, it can handle
tions for the potential.
even larger systems.
Athough ions are excluded from the interior of the
membrane by hard walls in LEMC, the electric field is
stillpresentthereandcouldbecomputed. Thedielectric
II. MODELS AND METHODS
constantisthesamethereasintheelectrolyte(ǫ=78.5),
therefore,the surfaceof the membraneis notadielectric
A. Models boundaryandpolarizationchargesarenotinducedthere
(Fig. 1A). In the case of PNP, the interior of the mem-
When we extract macroscopic information (currents, braneisnotpartofthecomputationaldomain(Fig.1B).
profiles, etc.) from a microscopic model, we construct An appropriate Neumann boundary condition is applied
a model that contains the interactions between the par- on the surface of the membrane in order to mimic the
ticles and the external constraints (hard walls, applied system used in LEMC. Boundary conditions as handled
field, etc.). This is equivalent with defining the Hamil- in the two methods are detailed in the following subsec-
tonian of the system precisely. This model then can tions.
be studied with different statistical mechanical methods
(simulations or theories). Whether a disagreement with
experimentsisduetooversimplificationsinthemodeling
B. Poisson-Nernst-Planck theory
or the approximations in the method can be sorted out
by comparing to particle simulations, where approxima-
tionsinthemethodareusuallyabsent(systemsizeerrors Introduced for modeling semiconductors,49,50, PNP
and statistical noises are still present). was soon adapted to the modeling of biological ion
Although the separation of model and method is not channels32,33,51–60 as well as synthetic nanopores61–63.
so distinct in PNP, we describe the two models together The classicalPNP is aself-consistentsystemprovidinga
in this subsection in order to emphasize similarities and flux that satisfies the continuity equation
differences. Note, however, that the term “method” in
the case of PNP refers to the physical equations used in ∇·ji(r)=0. (9)
PNP (see next subsection), not the numerical method
with which we solve the PNP equation. The flux is computed from the NP equation (Eq. 1) for
The electrolyte is modeled in the implicit solvent the linear and Eq. 6 for the nonlinear version) with the
framework in both cases. Water is a continuum back- electrochemicalpotential, µ ,definedinEqs.3and7,re-
i
ground,whose energeticeffect is takeninto accountby a spectively. Themeanelectrostaticpotential,Φ(r),iscon-
dielectric screening (ǫ = 78.5 in the denominator of the nectedtotheconcentrationprofilesinbothcasesthrough
4
(A) Γ Γ (B) Γ Γ
Φ=0 L 78.5 R Φ=U Φ=0 L =g R Φ=U
c=cL --ε=-+++ c=cR c=cL En ΓW c=cR
---+++
P N
5 σ σ Γ
ε=78.5 =78. ε=78.5 E=nE=n M
ε
FIG.1: Geometryofcomputationdomain(A)intheNP+LEMCsystemand(B)inthePNPsystem. (A)Boundaryconditions
fortheNP+LEMCsystemareprescribedforthetwohalf-cylinders(darkandlightgreenlines,ΓL andΓR domains)onthetwo
sidesofthemembrane. Dirichletboundaryconditionsareappliedbyusingtheappropriateappliedpotentialobtainedbysolving
theLaplace equation (alinearinterpolation way usedinsidethemembrane, seethedottedlines). Boundary conditionsforthe
concentrations are ensured by using appropriate electrochemical potentials at the boundaries that correspond to the chemical
potentials producing the prescribed concentrations. The domains outside the green lines are in thermodynamic equilibrium,
wherethechemicalpotentialisconstant,soequilibriumGCMCsimulationsareperformedthere. Porechargesarefreecharges
present explicitly in the simulation cell. They are placed on the pore wall on a grid as partial point charges. The dielectric
constantisthesameeverywhere,includingtheinteriorofthemembrane. (B)ThePNPcomputationalcellexcludestheinterior
of the membrane from the solution domain. The pore charges are polarization charges that are induced as a result of the
prescribed Neumann boundary conditions on the pore wall (red and blue lines, ΓW). On the surface of themembrane (brown
lines, ΓM), a Neumann boundary conditions is applied in order to mimic the NP+LEMC solution. On the two half cylinders,
thesame boundaryconditions are used as in NP+LEMC (ΓL and ΓR).
the Poisson equation: boundary conditions
e ci(r)=cLi and Φ(r)=0 on ΓL
−∇·(ε∇Φ)= z n (r). (10)
ǫ0 Xi i i ci(r)=cRi and Φ(r)=U on ΓR (11)
The third part are the regions of the membrane which
Here n (r) is the number density of ions (measured in
i areattachedtothe bathsandaredenotedbyΓ (brown
m−3)connectedtoconcentration(measuredinmol/dm3) M
lines in Fig. 1B). As the membrane is impenetrable for
through n (r)=1000N c (r). We further denote by N
i A i A the particle flux, we set the flux to be equal to 0 there.
Avogadro’snumberandbyǫ thepermittivityofvacuum.
0 InLEMCsimulationsthemembraneispenetrableforthe
Equation 3 is equivalent to the application of Boltz- electric field, whichis notthe case inPNP.Therefore we
mann’s distribution, which provides the statistical me- impose the boundary conditions
chanical description of the system. Together with the
Poisson equation it forms the PB theory. Coupling ∂Φ(r)
j (r)·n =0 and =g(r) on Γ , (12)
these to the NP equation applies the theory for a non- i M ∂n M
M
equlibriumsituationandprovidesthePNPtheory. From
a physical chemical point of view, Eq. 3 corresponds where n is the outer normal on Γ and the function
M M
to the statement that the electrolyte solution is ideal g(r) is supposed to mimic the LEMC case (where there
(µBiMF(r) = 0). While many theories were developed isanelectric fieldacrossthe membrane). Moreprecisely,
in the last years to overcome this limitation (e.g. DFT), itisobtainedbysolvingaLaplaceequationwithzeroleft
hereweemployavariantofPNPthatfeaturesnon-linear hand side without permanent charges and with bound-
cross-diffusion in the continuity equations which limits
ary condition Eq. 11 in the domain of Fig. 1A. Then,
the maximal density and thus takes into account the fi- evaluating the normal derivative of this solution at the
nite size of the ions. boundary Γ yields the function g(r). This additional
M
The (n)PNP systems are solved inside the computa- Neumann boundary condition matches the value of ap-
tional domain, whose boundary is separated into four plied potential crossing the membrane in the LEMC.
parts as shown in Fig. 1B. The first two parts corre- The last part of the boundary is on the inside wall
spond to the left and righthalf-cylinders (dark and light of the pore, called Γ . As it is a part of the membrane,
W
green lines in Fig. 1B) and are denoted by Γ and Γ . whichisimpenetrablefortheparticles,no-fluxconditions
L R
TheseregionsarethesameinNP+LEMC.Boththecon- arealsoimposedforthecurrent. Thepermanentcharges
centration and the potential are set using the following induce an additional electric field and are included by
5
another Neumann boundary condition: whereVk isthevolumeofBk,Nk isthenumberofparti-
i
clesofcomponentiinBk beforeinsertion,∆Uk istheen-
∂Φ(r) ergy change associated with the insertion (including the
j (r)·n =0 and =σ (z) on Γ , (13)
i W ∂nW 0 W effect of the externalfield), and µki is the configurational
(total minus µ0) electrochemical potential of component
i
where σ0 = σ and σ0 =−σ for z <0 and z >0, respec- iinBk. Inthe particledeletionstepwerandomlychoose
tively, and nW is the outer normal on ΓW. aparticleofcomponentiinsub-volumeBk anddeleteit.
Oneofthemostpopularsimplificationofthe fullPNP
The deletion is accepted with probability
model is the 1D reduction. It is broadly used and com-
pared with experimental data61. It also allows to sim-
Nk −∆Uk−µk
ulate very long nanopores with complex geometry using pk =min 1, i exp i . (16)
reasonablecomputationaltime. Thederivationofthe1D i,DEL (cid:26) Vk (cid:18) kT (cid:19)(cid:27)
model is based on the assumption that the length of the
nanopore is significantly larger than its radius. Since in Here, Nk is the number of particles of component i in
i
oursetup this is notthe case,weperformallsimulations subvolume Bk before deletion.
in two spatial dimensions. The energy change ∆Uk contains the effect of the full
To actually solve the 2D (n)PNP system we use the simulationdomainoutsidesubvolumeBkincludingshort-
well-known Scharfetter–Gummel scheme which is based range interactions such hard-sphere exclusions between
on a transformed formulation of the system in exponen- ions and hard-wall exclusion with membrane wall. The
tialvariables,see64 fordetail. Weusea2Dfiniteelement configurational space is sampled properly, because the
method for the actual implementation and a triangular ionsexperiencethepotentialproducedbybillionsofpos-
mesh containing 20−60 thousand elements, depending sible configurations, not just a mean potential as in the
on the radius of the pore. The mesh is also non-uniform case of PNP.
in order to obtain high accuracy, especially close to the
The effect of the applied potential is also included in
pore entrances. ∆Uk. The applied potential is computed by solving the
Laplace equation with the Dirichlet boundary condition
ofEq.11 forthe boundarysurfaceconfining the solution
C. Nernst-Plank equation coupled to Local
domain. Theboundaryconditionsforconcentrations(see
Equilibrium Monte Carlo
Eq.11)aresetbyfindingtheappropriatechemicalpoten-
tialsinthetwobathsthatproducethedesiredconcentra-
To solve the NP+LEMC system, an iterative proce-
tions in the GCMC simulations. We used the Adaptive
dure is needed, where µi is updated until the continuity GCMCmethod65 todeterminethesechemicalpotentials.
equation(Eq. 9) is satisfied. The procedure canbe sum-
The result of the simulation is the concentration ck
marized as i
in every volume element. The values ck and µk are as-
i i
LEMC NP ∇·j=0 signed to the centers of the volume elements and so the
µ [n] −−−−→ c [n] −−→ j [n] −−−−→ µ [n+1]. (14)
i i i i corresponding profiles are constructed. Both ck and µk
i i
fluctuate duringthe iterationprocess,sothefinalresults
The electrochemical potential for the next iteration,
are obtained as running averages.
µ [n+1],iscomputedfromthe resultsofthepreviousit-
i
eration, c [n], in awaythat they together producea flux The NP+LEMC technique has been applied to study
i
(through the NP equation) that satisfies the continuity transport through membranes1,2, calcium channels3,4,
equation. Details on the algorithm can be found in our and bipolar nanopores48.
original paper1.
The concentration profile in an iteration, c [n], corre-
i
spondingtotheelectrochemicalpotentialprofile,µ [n],is
i
III. RESULTS AND DISCUSSION
obtainedfromLEMCsimulations. Wedividethecompu-
tational domain (inside the green lines in Fig. 1A) into
volume elements and assume local equilibrium in these The reference point of all simulations corresponds to
volume elements. We assume that these local equilibria the following parameter set: voltages ±200 mV (200
can be characterized by local electrochemical potential mV is the ON, while -200 mV is the OFF state of the
values. We also assume that the gradient of the µ pro- nanopore), concentrations c = 0.1 and 1 M, surface
i
file defined this way is the driving force of ion transport charge σ = 1 e/nm2, and nanopore radius R = 1 nm.
as described by the NP equation (Eq. 1). Then, we vary the parameters systematically by chang-
TheheartoftheLEMCsimulationisaMCstep,where ing only one and keeping the others fixed. Rectification
weinsert/removeanioninto/froma volumeelementBk. is defined by |I(U)/I(−U)|, i.e. the ratio between the
The acceptance probability of an insertion is currents in the ON and OFF state, respectively. In our
case, this implies that it is always larger than 1. In all
Vk −∆Uk+µk figures we plot the NP+LEMC, PNP, and nPNP results
pk =min 1, exp i , (15)
i,INS (cid:26) Nk+1 (cid:18) kT (cid:19)(cid:27) with symbols, solid lines, and dashed lines, respectively.
i
6
R = 1 nm, σ = 1 e/nm2 c = 1 M, R = 1 nm, U = ±200 mV
NP+LEMC 200 NP+LEMC
40 Rectification 100 PNP ON
lin. PNP
nPNP
non lin. PNP
30
50
A 150
p
I / 1200 0 c = 100.01 M 2000 | I / pA | 100 ectification112050
0 R 5
50 0
200 0 0.5 1 1.5
Rectification c = 1 M
10 OFF
150
0
5 0 0.5 1 1.5
A 100
p σ / e/nm2
I / 50 0 100 2000
FIG. 3: The absolute value of the current as a function of
0
σ (characterizing the strength of the polarity of the pore) in
the ON and OFF states (200 vs. -200 mV, respectively) as
-200 -100 0 100 200 obtained from NP+LEMC, PNP, and nPNP (symbols, solid
U / mV curves, dashed curves, respectively). The inset shows rectifi-
cation. The model parameters are c=1 M and R=1 nm.
FIG. 2: Current-voltage curves for concentrations c=0.1 M
(top panel) and c = 1 M (bottom panel) as obtained from voltages of opposite signs are the same and rectification
NP+LEMC, PNP, and nPNP (symbols, solid curves, dashed is 1. Figure 3 shows current values in the ON and OFF
curves, respectively). The insets show rectification as com-
states as functions of σ. As σ is increased, the current
puted from theratio of the ON and OFF state currents (the
increases in the ON state, while decreases in the OFF
absolute values). The model parameters are R = 1 nm and
σ=1 e/nm2. state. Rectification, therefore, improves as the strength
of the polarity of the pore increases. The σ-dependence
is well described by (n)PNP qualitatively. The errors
manifest in the fact that rectification is underestimated
A. Comparison of I-U curves and rectification
by (n)PNP.
behavior
One source of the errors is that the effective cross sec-
tion of the pore through which the centers of ions can
First,welookatthenanoporeasadevicethatgivesan moveissmallerinthecaseofthechargedhardsphereions
output signal (current) as an answer to the input signal used in LEMC (R−0.15 nm, where 0.15 nm is the ionic
(voltage). Therelationoftheseisthetransferfunctionof radius) than in the case of point ions used in (n)PNP
the device. Then, we study various profiles (concentra- (the whole pore radius, R, is used in (n)PNP). (n)PNP,
tion,potential,chemicalpotential)andtrytounderstand therefore, systematically overestimates current in both
the differences between PNP and NP+LEMC. the ON and OFF states as seen in Fig. 3. The overesti-
Figure2showscurrent-voltage(I-U)curvesforthecon- mationofthedenominator(OFFcurrent)dominatesthe
centrationsc=0.1and1M.Rectificationisobservedus- ratio. Rectification, therefore, is underestimated.
ingallthethreemethods: thecurrentislargeratpositive One way to partially overcomethis difference between
voltages than at negative voltages (note that electrical the twomodelswouldbe usingthe effectivecrosssection
currents are multiplied with -1 in order to get positive of the finite ions (R−0.15nm) in the PNP calculations.
currents for positive voltages). Rectification increases In this case, Fig. 3 would show better agreement, but
with increasing |U| as shown in the insets. Agreement cause other problems, such as the presence of ions with
between NP+LEMC and (n)PNP data is better at low different diameters. Therefore, we decided to keep the
concentration (0.1 M) and smaller voltages as expected. pore cross section in PNP in this study as it is (R), but
The data fromnPNP areslightly better than those from point out the problems with this approach.
PNP, especially for c=1 M. Figure 4 shows the currents as functions of the elec-
The value of the σ parameter can be considered as a trolyteconcentration,c. Currentsdecreasewithdecreas-
measureofthe nanopore’spolarity. Atσ =0e/nm2,the ing c as expected, but the current decreases faster in
pore is unchargedand symmetric, so currents at the two the OFF state, so rectification increases with decreas-
7
R = 1 nm, U = ±200 mV, σ = 1 e/nm2 c = 1 M, U = ±200 mV, σ = 1 e/nm2
1500 NP+LEMC
500
40 PNP
220500 234000000Rectification ON 1000 2300Rectificaion nPNP
A | 100 A | 10
p 150 0 p ON
| I / 0.01 0.1 1 | I / 1 2 3
100 NPNP+PLEMC OFF 500
nPNP
50
OFF
0 0
0.01 0.1 1 1 2 3
c / M R / nm
FIG. 4: The absolute value of the current as a function of FIG. 5: The absolute value of the current as a function of
theelectrolyteconcentration in theONand OFFstates(200 theporeradius intheONandOFFstates (200 vs.-200mV,
vs. −200 mV, respectively) as obtained from NP+LEMC, respectively) as obtained from NP+LEMC, PNP, and nPNP
PNP, and nPNP (symbols, solid curves, dashed curves, re- (symbols,solidcurves,dashedcurves,respectively). Theinset
spectively). Theinsetshowsrectification. Themodelparam- showsrectification. ThemodelparametersareR=1nmand
eters are c=1 M and R=1 nm. σ=1 e/nm2.
B. Analysis of profiles for concentration, electrical
potential, and electrochemical potential
ing concentration, a well-known result. The explanation
Togetadditionalinsightsintothephysicalmechanisms
is that depletion zones dominate the currents in bipolar
beyondthedevice-levelbehavior,wealsoanalyzeprofiles
nanopores,but depletion zones are more depleted at low
for the concentration, electrical potential, and electro-
concentrations. Changing the sign of the voltage from
chemical potential.
positive (ON) to negative (OFF), therefore, can deplete
In Fig. 6, we plot the concentration profiles for c = 1
thedepletionzonefurthermoreefficientlyatlowconcen-
(panel A) and 0.1 M (panel B) in order to study the
trations.
differences between high and low concentrations. This
figureshowsthe resultsforσ =1e/nm2. Figure7shows
AgreementbetweenNP+LEMCand(n)PNP is better the same concentration profiles but for σ =0.25 e/nm2.
in the ON state. The nonlinear version of PNP works The curves show that the ions have depletion zones in
better in this case, because it handles crowding better. the middle of the pore and in the zone, where they are
In the OFF state, (n)PNP systematically overestimates theco-ions(havingionicchargewiththesamesignasthe
thecurrentpartlyfromthereasondiscussedabove. Rec- pore charge, σ). We distinguish basically four regions:
tification, interestingly, is underestimated by (n)PNP at
1. left bath, near the membrane (z <−3 nm)
large, while overestimated at small concentrations.
2. theleftpartoftheporewithpositivesurfacecharge
Finally, we show the dependence of currents on the (−3 < z < 0 nm, N region) – anions the counter-
pore radius in Fig. 5. Currents increase with widen- ions, cations the co-ions
ing pores as expected. The relative difference between
3. the right part of the pore with negative surface
the ON and OFF states decreases as R increases. Rec-
charge (0 < z < 3 nm, P region) – cations the
tification is the result of the interplay between the ef-
counter-ions, anions the co-ions
fect of pore charges and the applied potential. The av-
erage distance of pore charges from the ions increases
4. right bath, near the membrane (z >3 nm)
as R increases, therefore, the pore charges get less and
less able to produce the depletion zones inside the pore. Intheaccessregions,closetotheporeentrances(regions
(n)PNP qualitatively reproduces the behavior obtained 1 and 4) ionic double layers are formed. Double layer is
from NP+LEMC. Also, the systematic underestimation common name for the separation of cations and anions
of rectification is present for all pore radii studied. (polarization of the ionic distributions) as a response to
8
(A) σ = 1 e/nm2, c = 1 M (B) σ = 1 e/nm2, c = 0.1 M
5
Symbols: LEMC ON Symbols: LEMC ON 3
Solid: PNP Solid: PNP
4 Dashed: nPNP Dashed: nPNP
M anion cation 2
) / 3
z
(
c 2 anion cation
1
1
0
2 anion cation 1
anion cation
M
0.1
) /
z
( 1
c
0.01
OFF OFF
0.001
0
-4 -2 0 2 4 0
z / nm z / nm
FIG. 6: Concentration profiles of cations and anions as obtained from NP+LEMC, PNP, and nPNP for (A) c=1 M and (B)
c = 0.1 M for parameters R = 1 and σ = 1 e/nm2. These concentration profiles have been computed by taking the average
number of ions in a slab and dividing by the available volume. For −3<z <3 nm, the cross section of the pore was used to
obtain this volumein both methods.
thepresenceofachargedorpolarizedobject. Inthiscase, brane. Atleast,thisseemstobesuggestedbytheresults
double layers appear partly as a response to the applied of NP+LEMC.
field, partly as a responseto the chargeimbalanceinside The double layers have opposite signs in the ON and
thepore. Realizethatthesignofthedoublelayer(which the OFF states that can be explained through the mean
ionsaretheco-ionsandcounter-ionsinthe doublelayer) electrical potential profiles that have two components
depends on the sign of the applied voltage. produced by all the free charges,ΦFREE(r), and induced
The basic reason of rectification is that the ions are charges, ΦAPP(r), in the system. In this study, induced
more depleted in their depletion zones in the OFF state charges appear at the boundaries where the boundary
thaninthe ONstate;cationconcentrationintheNzone conditions are applied, therefore, they produce the ap-
islowerintheOFFstatethanintheONstate,forexam- plied potential, ΦAPP(r). The total mean potential,
ple. Basically,thedepletionzonesarecausedbythepore therefore, is obtained as
charges. The applied field modulates the effect of pore
Φ(r)=ΦFREE(r)+ΦAPP(r). (17)
charges, therefore, it increases or decreases concentra-
tionscomparedtothezero-voltagecase. Depletionzones InthecaseofNP+LEMC,thedoublelayersarenecessary
are the main determinants of the current, because they toproducetheΦFREE(z)componentthatcounteractsthe
arethehigh-resistanceelementsofthesystemmodeledas appliedfield,ΦAPP(z). Figure8Ashowsthattheslopeof
resistorsconnectedinseriesalongthe ionicpathway. So, ΦFREE(z) is the opposite to the slope of ΦAPP(z) in the
if depletion zones are more depleted, current is reduced. bulks,sotheirsum(TOTAL)hasthe slopeclosetozero.
It is important,however,thatnot onlythe co-ioncon- This is necessary because the bulks are low-resistance
centrations decrease by switching from ON to OFF, but elements, where the potential drop is small.
also the counter-ionconcentrations. As a matter of fact, In the case of (n)PNP, this phenomenon depends on
this is crucial, because co-ions are brought into their de- the imposed boundary conditions, Eq. 12, on the mem-
pletionzoneswiththehelpoftheirstrongcorrelationsto brane surface. Using, for example, g = 0 yields to-
counter-ions. So there are less co-ions because there are tally different results which are in poor agreement with
less counter-ions. The quantity of counter-ions, on the NP+LEMCasfaras the structureofthese double layers
other hand, seems to be related to the double layers at is concerned (the behavior inside the pore is less influ-
the entrances of the pore on the two sides of the mem- enced).
9
(A) σ = 0.25 e/nm2, c = 1 M (B) σ = 0.25 e/nm2, c = 0.1 M
0.8
ON ON
anion cation
anion cation
0.6
1.5
M
) / 0.4
z
(
c 1
Symbols: LEMC
0.2
Solid: PNP Symbols: LEMC
Dashed: nPNP Solid: PNP
0.5
0
1.2
anion cation anion Dashed: nPNP cation
0.2
1
M 0.15
0.8
) /
z
c( 0.6 0.1
0.4 0.05
OFF OFF
0.2 0
-4 -2 0 2 4 -4 -2 0 2 4
z / nm z / nm
FIG. 7: Concentration profiles of cations and anions as obtained from NP+LEMC, PNP, and nPNP for (A) c=1 M and (B)
c=0.1 M for parameters R=1 and σ=0.25 e/nm2.
Comparingthecounter-ionprofilesinthedoublelayers and7). In the OFF state, (n)PNP usually overestimates
and in the neighbouring half nanopores (Figs. 6 and 7), concentrations causing the overestimation of current as
wecanseethatiftherearelesscounter-ionsinthedouble we have seen before. This is counterintuitive, because
layer,therearelesscounter-ionsinthehalfnanoporetoo it was said that (n)PNP is better at low concentrations,
(see anions on the left hand side in the OFF state com- but pore concentrations are higher in the ON state. We
pared to the ON state, for example). Although the de- canresolvethiscontradictionifweconsiderthatthesys-
crease of counter-ionconcentrationin the pore is related tem’sbehaviorisaresultofthebalanceofbasicallythree
tothedecreaseoftheconcentrationofthesameioninthe effects: (1) interaction with the fixed pore charges, (2)
neighbouring double layer,it would be an overstatement interaction with the fixed applied field, and (3) mutual
to say that one is a consequence of the other. and complicated interactions between ions. The mutual
weight of these terms is different in the ON and OFF
Rectification works without this coupling between ion
states.
quantities in the double layer and in the nanopore. For
example, rectification is reproduced in the case of PNP In the ON state, pore charges and applied field act
with boundary condition g = 0 although with worse in the same direction, so they dominate the energy and
agreementwithNP+LEMC.Furthermore,theformation errors in the ion-ion term have less effect. In the OFF
ofthedoublelayersisabsentinMDsimulationsusingex- state, however, pore charges and applied field act in the
plicit water, still, rectification is present. MD results us- opposite directions, so their sum is smaller and the ion-
ingexplicitwaterareingoodagreementwithNP+LEMC ion term has a larger weight and the BMF term with
results using implicit water48. These contradictions re- it.
quire more study, but it seems that the formation of the
Ournextgoalistobetterunderstandthedifferentcon-
double layers is rather related to boundary conditions
tributionsofthecomponentsofthe totalelectrochemical
and larger-scale effects, while the structure of the ionic
potential µ as defined in Eqs. 2–5. The electrical com-
i
profiles inside the pore is rather related to local effects ponent µEL is defined as the interaction with the (to-
i
such as interaction with pore charges, applied field, and
tal) mean electrical potential that is shown in Fig. 8A.
other ions.
Note that the BMF term is fully included in the CH
As far as the agreement between the NP+LEMC and termandtherefore,inthe caseof(n)PNP,µCH(r)is just
i
the theoretical profiles is concerned, it is generally bet- µID(r), while it also contains the BMF term in the case
i
ter in the ON state than in the OFF state (see Figs. 6 of NP+LEMC.
10
(A) σ = 1 e/nm2, c = 1 M (B) σ = 1 e/nm2, c = 1 M
0.2 ON TOTAL 2
APP T
0.15 Symbol: NP+LEMC 0 n / k
o
V Lines: PNP -2 ni
mponents / 0.00.51 FREE µµµµICBEDMLH--µFµCCHH((LL)) ON --64 onents of a
co 0 µ−µCH(L) -8 mp
al co
nti al
al pote 0.005 FREE µµµ−CELHµ -C(µPHCN(HLP()L )()P (NPPN)P) 68 potenti
ectric-0.05 OFF 4 mical
El e
-0.1 2 ch
-0.15 APP 0 ctro
e
OFF El
-0.2 TOTAL -2
-6 -4 -2 0 2 4 6 -4 -2 0 2 4
z / nm z / nm
FIG. 8: (A) Electrical potential profiles and components (see Eq. 17) as obtained from NP+LEMC and PNP. Component
ΦFREE(z) is the product of ions and pore charges in the system, while ΦAPP(z) is the applied potential computed from the
Laplace equation with Dirichlet boundary conditions. (B) Electrochemical potential profiles and components (see Eqs. 2, 3-5)
as obtained from NP+LEMC and PNP. The ideal (µIiD(z) = µ0i +kTlnci(z)), the electrochemical (µi(z)), and the chemical
(µCiH(z)) terms are shifted to zero by deducting µCiH(L), which is the value of the chemical term in the left bath. In the case
of PNP the ID and CH terms are the same, so µBiMF =0. Results are shown for the anion; data for the cation do not reveal
new insights (subscript iis dropped in the legend). Parameters are c=1 M, σ=1 e/nm2, and R=1 nm.
Figure8Bshowsthefullelectrochemicalpotentials,the hand, indicates that there is an error in “chemistry”, so
CHterms,andthe ELterms. Inthe caseofNP+LEMC, there is a possibility for further errors in both the c (r)
i
we also plot the lnc (z) term (denoted as ID) and the and Φ(r) profiles inside the pore. Those potential errors
i
BMF term. The ID and CH terms, as well as the total can eventuate inside the pore and become visible in all
electrochemical potential, are all shifted by the value of profiles. Local fluctuations in the BMF term inside the
the CH term in the left bath (µCH(L)). In this way, the pore indicate how seriously do the errors of the mean-
i
µ (z), µCH(z), and µEL(z) contributions take the value field treatment of PNP contribute to inaccuracies of all
i i i
zero at the left edge of the plot. the profiles inside the pore.
The errors in µ have three components: the error in
i
reproducing (1) the lnc term, (2) the EL term, and (3)
i
the BMF term that can be identified with errors in re- IV. SUMMARY
producing the particle correlations which are missing in
PNP,duetothemeanfieldapproximation. Thefirsttwo
The general conclusion is that the BMF term is small
errorshavedifferentsignsandtendtobalanceeachother.
andtheagreementbetweenPNPandNP+LEMCisquite
TheyarecoupledthroughthePoissonequation,sointhe
good. Yet, since the mean field theory does not capture
limiting case of agreeing c profiles, the Φ profiles agree
i the OFF state behavior as good as the ON state be-
if the boundary conditions are also the same.
havior, derived quantifies as the rectification cannot be
Inthis case,the NPequationwouldgivethe sameflux predicted that well. Still, the results are very promising
if the BMF term were constant, because ∇µ would be given that these calculations have been performed for a
i
the same in the two methods. Therefore, the real source narrow (R = 1 nm) and short (6 nm) pore with experi-
of errorsis not the magnitude of the BMF term, but the mentally typical, but quite large surface charges (σ ∼ 1
r-dependence of the BMF term, that is, the fact that e/nm2). This indicates that the 2D PNP used in this
ionic correlations are different inside the pore than out- study is an appropriate tool to study more realistic ge-
side. The nonzero value of the BMF term, on the other ometries (wider and longer pores), at least, as far as the
