Multiscale 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