Simulating Ion Transport with the NP+LEMC Method. Applications to Ion Channels and Nanopores.

Abstract We describe a hybrid simulation technique that uses the Nernst-Planck (NP) transport equation to compute steady-state ionic flux in a non-equilibrium system and uses the Local Equilibrium Monte Carlo (LEMC) simulation technique to establish the statistical mechanical relation between the two crucial functions present in the NP equation: the concentration and the electrochemical potential profiles (Boda, D., Gillespie, D., J. Chem. Theor. Comput., 2012 8(3), 824–829). The LEMC method is an adaptation of the Grand Canonical Monte Carlo method to a non-equilibrium situation. We apply the resulting NP+LEMC method to ionic systems, where two reservoirs of electrolytes are separated by a membrane that allows the diffusion of ions through a nanopore. The nanopore can be natural (as the calcium selective Ryanodine Receptor ion channel) or synthetic (as a rectifying bipolar nanopore). We show results for these two systems and demonstrate the effectiveness of the NP+LEMC technique.


Introduction
We dedicate this paper to honoring Professor János Liszi, who was the supervisor of one of the authors (DB) and respected senior to another (MV), supporting to the careers of both. We (DB and VM) publish this paper together with our students (DF and EM) to demonstrate that the lesson of Professor Liszi -one of the most important goals of senior scientists is to pave the way to the junior ones -has not been left unconsidered.
The goal of this work is to present a computer simulation technique for computing steady-state ion transport that was developed and applied in our research group with the essential help of Dirk Gillespie [1]. The method is based on the Nernst-Planck (NP) transport equation coupled to the Local Equilibrium Monte Carlo (LEMC) simulation technique and is described in Sec. (2) in detail.
The method, called NP+LEMC, was applied for various problems in the last couple of years. These problems include particle transport through model membranes [1,2], ion channels [3][4][5], and nanopores [6][7][8]. The basic goal is to have a computationally efficient simulation technique using reduced models of steady-state systems, where particles (ions being the focus) diffuse as a result of a maintained driving force that is a concentration gradient and/or an electrical potential gradient (that is, an electrochemical potential gradient).
We tend to regard these kind of systems as nanode- * Correspondence: fertig.david92@gmail.com vices that yield a macroscopic output signal (electrical current, for example) as a response to an input signal (voltage, for example). Reduced models have the advantage of including those degrees of freedom of the manyparticle system that are essential in reproducing device behavior, namely, the relationship of the output and input signals [6].
We present two case studies here to demonstrate both the effectiveness (and also possible source of errors) of the method and the power of reduced models. In the Ryanodine Receptor (RyR) calcium release ion channel, we needed to model only those amino acids of the channel protein that are inside or near the selectivity filter and hang into the ionic pathway. In the bipolar nanopore, we needed to reproduce the variation of the electric field along the pore axis properly. The most important reduction in the degrees of freedom, however, is that we model water as a dielectric continuum. The effect of using implicit solvent instead of an explicit solvent model was discussed in our previous work [6] through comparisons to molecular dynamics (MD) simulations.
Our computational method is not the only one that is able to determine ionic current through biological and synthetic nanopores using reduced models. While the Brownian Dynamics (BD) simulation method [9][10][11][12][13][14] is the obvious candidate to simulate this problem, it has disadvantages from the point of view of sampling the flux of ions. The Poisson-Nernst-Planck (PNP) theory [11,[15][16][17][18] also uses the NP equation to compute flux, but uses a mean-field approach, the Poisson-Boltzmann (PB) theory, to describe the statistical mechanical relationship between the concentration and electrochemical potential profiles.
To use something more state-of-the-art to provide this relationship beyond the PB level is the essence of our approach. We suggest the LEMC technique, which is a particle simulation method based on a three-dimensional model producing all the ionic correlations missing from PB. Gillespie et al. proposed a different technique, in which a density functional theory (DFT) was coupled to the NP equation [19,20] and used to method (NP+DFT) to describe the behavior of the same ion channel that we study here, the RyR ion channel [21][22][23]. DFT is a continuum theory, but a very sophisticated one, where ionic correlations are included through a series expansion of the free energy functional and the finite size of ions is taken into account through the fundamental measure theory [24]. Its accuracy was demonstrated by computing electrical double layers in extreme conditions and comparing to Monte Carlo (MC) simulations [25,26]. The disadvantage of this method is that it can be used efficiently only in one dimension. This reduction in dimensionality, however, is often feasible provided that we reduce the model in an intelligent way.

Modeling and methodology
A good model should be able to explain complex phenomena in a simplified way while it stays in agreement with experimental data. The nanopore systems (both natural and synthetic) studied in this work have some common features. The primitive model of electrolytes, that represents ions as charged hard spheres, is used. The interaction potential between two ions is defined by Coulomb's law in a dielectric material: where R i and R j are the radii, q i and q j are the charges of the different ionic species, 0 is the permittivity of vacuum, = 78.5 is the relative permittivity of the solvent, and r is the distance between two ions. Solvent (water, in this work) is not modeled explicitly; it is rather represented by two response functions. Energetically, solvent acts as a dielectric background that screens the electric field of ions by just dividing by in the Coulombpotential (Eq.(1)). Dynamically, solvent molecules collide with ions and impede their diffusion; we take that effect into account by a diffusion coefficient function, The ions of the electrolyte diffuse through a pore between two bulk containers (see Fig.1). The pore penetrates a membrane that separates the two containers. The pore and the membrane are also modeled in a reduced way; they are confined by hard walls with which the three-dimensional model is obtained by rotating the figure about the z-axis. The domain confined by the blue line is the non-equilibrium transport region (the solution domain of the NP+LEMC system) that includes the pore and the two access regions. The electrochemical potential in this region is not constant, but changes from one constant value to another (the values in the two bulk regions) in a monotonic manner to provide the driving force of ionic transport. Parameters R and H characterize system size.
the hard sphere ion cannot overlap. Other details of the pores (amino acid side chains and surface charges) will be presented in Sec.(3) for the respective ion channel and nanopore systems. Ion transport is steady-state; concentrations and electrical potentials, therefore, are kept fixed at the boundaries of the two baths on the two sides of the membrane (blue line in Fig.1). We assume that ion transport is described by driftdiffusion. Our procedure is a hybrid method that separates the configuration degrees of freedom of ions (ion positions) from the kinetic degrees of freedom (ion velocities). In BD simulations they are treated together and ionic current is measured by counting ions that pass through the pore. The BD method has the disadvantage of weak sampling of flux especially at low concentrations (as in our ion channel example) and in cases, when current is limited by ionic depletion zones (as in our bipolar nanopore example).
Separation of the two parts of the phase space is a usual practice in equilibrium statistical mechanics, where the kinetic part is described by the ideal gas equations, while the excess quantities that produce all the peculiarities beyond ideal-gas behavior are from the potential energy, that, in turn, depends only on the configuration coordinates [27][28][29][30]. Out of equilibrium, however, this separation is not so obvious. Classical mechanics is still valid, so particles' motion can be described by Newton's equa-tions of motion (MD simulations can be used). In the continuum solvent framework, Langevin's equation is used in BD simulations.
The meaning of various thermodynamic quantities, however, loses the solid ground that it has in equilibrium. There is no well-established non-equilibrium statistical mechanics finding its way into textbooks despite the fact that many authors made considerable advances in this field [31][32][33][34][35][36]. The most problematic quantities are those containing entropic effects, such as the chemical potential, µ i . In the case of charged particles, we use the term electrochemical potential, but we talk about the same thing. This is a crucial quantity, because its homogeneity defines thermodynamic equilibrium (let us assume that the other two intensive parameters, temperature and pressure, are constant). If µ i is not constant, species i will diffuse until it becomes constant according to the second law of thermodynamics. The gradient of µ i (r), therefore, is the driving force of particle diffusion.
The empirical transport equation that is most widely used to describe transport of charged particles (electrodiffusion) is the NP equation, where D i (r) is the diffusion coefficient profile of ionic species i, c i (r) is the concentration profile, µ i (r) is the electrochemical potential profile, j i (r) is the particle flux density, k is Boltzmann's constant and T = 298.15K is the temperature. The resulting j i (r) must satisfy the continuity equation: The NP equation separates the configuration and kinetic parts of the phase space in a way that it provides flux as a function of three profiles that, in turn, depend only on configuration coordinates. Thus, we reduced our statistical mechanical task to averaging over states in the configuration space just as we did in the case of equilibrium statistical mechanics. If we treat D i (r) as a parameter, the task is reduced to finding the proper relation between c i (r) and µ i (r). Defining local concentration in a nonequilibrium situation is the same as in equilibrium: we compute the average number of particles in a small volume element and divide it with the volume of that element (though we will also compute concentration in a different way on the basis of the Potential Distribution Theorem [37], see later). Eletrochemical potential, on the other hand, is more problematic. The trick in this case is that we assume local equilibrium (LE). If we divide the simulation cell into subvolumes B α , we can characterize this subvolume with the electrochemical potential, µ α i , and assume that this value is constant in B α . We designate the µ α i and c α i values to the mass centers of the B α volume elements. The assumption of LE is also present in other methods (though not necessarily stated explicitly), where the c i (r) vs. µ i (r) relationship is considered, such as in the PNP theory or in the NP+DFT method of Gillespie et al. [19,20].
The main proposal of our technique was to use an MC method for this purpose. We assumed that the subvolumes are open systems and they are in LE with a constant volume (V α ), temperature (T ), and electrochemical potential (µ α i ). Thus, we suggested to use Grand Canonical Monte Carlo (GCMC) simulations, where particle insertion/deletions are attempted in the simulation cell. The acceptance probability depends on which subvolume we try to insert to (or where is the particle that we try to delete): Here, N α i is the number of ions of type i in subvolume B α before insertion/deletion, ∆U (r) is the change of the system's potential energy during particle insertion to position r (or deletion from there), χ = 1 for insertion, and χ = −1 for deletion. The difference between this method and equilibrium GCMC is that µ i is space-dependent and the acceptance criterion is referred to a given subvolume instead of the whole simulation cell. The effect of the surrounding of subvolume B α , however, is taken into account in the simulation through the energy change that includes all the interactions from other subvolumes, not only interactions between ions in B α .
Although it is tempting to view the subvolume B α as a distinct thermodynamic system with its own ensemble of states and to include the effect of other subvolumes as an external constraint, this is not the case. The ensemble of states belongs to the whole system, because ion configurations in subvolume B α should be collected for every possible ion configurations of all the other subvolumes. Therefore, the independent variables of this ensemble are T and {V α , µ α i }, where α and i run over the volume elements and particle species, respectively. For comparison, the variables in global equilibrium are T , V , and µ i , where V is the total volume and µ i does not depend on space.
We solve the NP+LEMC system in an iterative way. The electrochemical potential is adjusted until conservation of mass (∇ · j i (r) = 0) is satisfied. The procedure can be summarized as , on the basis of the divergence-theorem (also known as Gauss-Ostrogradsky's theorem). The continuity equation is converted to a surface integral: where volume B α is bounded by surface S α and n(r) denotes the normal vector pointing outward at position r of the surface. Every S α surface is divided into S αβ elements. Along these elements, B α and B β are adjacent cells. It is assumed that the concentration, the gradient of the electrochemical potential, the flux density and the diffusion coefficient are constant on a surface S αβ . They are denoted by hat:ĉ αβ i , ∇μ αβ i ,ĵ αβ i , andD αβ i . Theĉ αβ i values are obtained from the values c α i and c β i via linear interpolation. The ∇μ αβ i values are also obtained from µ α i and µ β i assuming linearity. Thus the integral in Eq.(6) for a given surface S α is replaced by a sum over the surface elements that constitute S α : where a αβ is the area of surface element S αβ and n αβ is the outward normal vector in the center of S αβ . The iteration procedure is described by the following steps: 1. 4. To achive faster and more robust convergence in the case of large driving forces, the electrochemical potential used in the (n + 1)th iteration is mixed from the values calculated in the (n + 1)th iteration from Eq. (8) and the values mixed in the nth iteration: where b i is a mixing parameter that determines the ratio of mixing. If the parameter is close to 1, faster iteration can be achieved, however, it may result in the system fluctuating between local minima. Smaller b i values prevent this fluctuation at the price of making the convergence slower.
5. The input of the (n + 1)th LEMC simulation is The system of linear equations contains the boundary conditions for the electrochemical potential and concentration via the boundary elements that are at the system's boundaries (blue line in Fig.1). If the S αβ face is on the system's boundary, the values ofĉ αβ i andμ αβ i are those prescribed on that face. The electrochemical potentials on the left (L) and right (R) boundaries are computed from and where µ EX,L i and µ EX,R i are the excess chemical potentials in the absence of an external field (determined by the Adaptive GCMC method of Malasics et al. [38]), while q i Φ L and q i Φ R are the interactions with the applied electrical potentials. Prescribing Φ L and Φ R on the system's boundary means that we use an electrostatic Dirichlet boundary condition. Voltage is defined as U = Φ L − Φ R . Ultimately, we have boundary conditions for ion concentrations on the two sides of the membrane and for the voltage (ground is on the right in this study).
The energy change ∆U contains not only the interactions between particles (and interactions of particles with the pore), but also the interaction with an external electrical potential, Φ appl (r). This applied potential is calculated by solving Laplace's equation for the system (inside the blue line) with the prescribed Dirichlet boundary condition on the system's boundaries (Φ L and Φ R on the left and right blue line, respectively).
In the portion of the blue line inside the membrane, we interpolate between Φ L and Φ R . We solve this equation with the Induced Charge Computation method [39,40].
In the present geometry, the elementary cells are ∆z× ∆r rectangles in the (z, r) plane. In three-dimensional space, these correspond to concentric rings with the zaxis in their centers.
The concentration in an elementary cell B α can be computed from dividing the average number of ions in the cell with the volume of the cell. This route is disadvantageous if ion concentrations are small. An alternative method is based on the Potential Distribution Theorem [37] that practically corresponds to the Widom particle insertion method [41,42]. The total excess chemical potential (containing also the interaction with the applied field) is µ α i − kT ln c α i and can be computed as where ∆U α i is the energy change associated with the insertion of a test particle of species i into a randomly chosen position in the elementary cell α. This is the same energy computed in the particle insertion step of the LEMC technique, therefore, it does not require additional computational cost. The concentration can be calculated from Eq.(13) because µ α i is known (that is the input of the simulation) and the ensemble average on the right hand side can be obtained from the LEMC insertions. This route provides a more accurate value for c α i , because this sampling is not discrete as just counting particles, but rather continuous, because we always gain information from the simulation at every ion insertion via the energy ∆U α i . The route of counting particles, however, might be better, when concentrations are very large in the volume element. In the case of long enough simulations, the two methods give identical results (within a statistical error).
Because of this statistical error, the iteration does not converge to an exact value of the ionic flux density, but it fluctuates around a limiting value. The final solution is obtained from a running average over iterations. Longer LEMC simulations and more iterations result in a more reliable outcome. The net flux of the diffusing ions through the cross section of the pore is calculated via averaging: where R(z) denotes the radius of the pore at coordinate z (|z| < H memb /2, where H memb is the thickness of the membrane) and n z is the unit vector along the z-axis. The electrical current of an ion is then and the total current is

Calcium release channel
There are two important classes of calcium channels that have been our focus. The L-type calcium channel is found in excitable cell membranes [43] such as those of nerve cells and muscle cells (both cardiac and skeletal). These are strongly Ca 2+ selective channels, meaning that the presence of only a µM Ca 2+ in the bulk decreases the Na + current to half of its value compared to its value in the absence of Ca 2+ . This is called micromolar Ca 2+ affinity (or Ca 2+ block) because a very small amount of Ca 2+ is enough to block the channel [44]. The price of high Ca 2+ selectivity is low Ca 2+ current through this channel. The Ca 2+ ions flow into the muscle cell when the L-type calcium channels open in response to an electrical signal (action potential) in a small quantity that is not sufficient to initiate muscle contraction. The solution of Nature for this problem is an amplification mechanism.
A large amount of Ca 2+ ions are stored in organelles, the sarcoplasmic reticulum (SR), that are situated inside the muscle cells. There are calcium release channels (also known as RyR calcium channels) embedded in the membrane of the SR that are opened by the Ca 2+ ions provided by the L-type channel (in cardiac muscle) or via a direct link between the two types of calcium channels (in skeletal muscle). The RyR channels provide the large Ca 2+ flux that is necessary for muscle contraction by binding to tropomyosins on the actin filaments and allowing myosin to climb up the filament.
These channels are wider (also, less Ca 2+ selective) than the L-type channels. A large amount of experimental data is available for this channel [45][46][47]. On the basis of this database, Gillespie et al. developed a model for the RyR channel [21,22]. Although the model was a onedimensional reduction, studied with a one-dimensional NP+DFT, they were able to reproduce hundreds of experimental current-voltage curves.
Later, when the NP+LEMC method became available, we developed [4] a three-dimensional version of the model of Gillespie et al. The model is similar to that used by Boda et al. [3,4,40,[48][49][50][51][52][53][54] for the L-type calcium channel inspired by the "space-charge competition mechanism" proposed by Nonner et al. [55]. From the point of view of ion selectivity and permeation, the important part of the channel is the selectivity filter, the bottleneck of the pore. This region is lined by four P-loops that are parts of the four membrane-spanning subunits. These loops contribute four glutamic acids (in the L-type calcium channel) to the filter-lining region (the EEEE locus). The COO − groups of the carboxyl side chains are thought to reach into the pore lumen and interact with the passing ions [44]. In the case of the RyR channel, they are four aspartates (D4899, see Fig.2). Gillespie et al. [21-23, 46, 47] identified additional E and D amino acids in the two vestibules on the cytosolic and luminal sides (Fig.2).
It has been shown [3,4,40,[48][49][50][51][52][53][54][55] that the strong Ca 2+ selectivity can be reproduced by a reduced model where eight half-charged oxygen atoms (O 1/2− ) of the four COO − groups force a strict competition between Ca 2+ and Na + ions for space in the crowded filter. The exact positions of the O 1/2− ions are not known and even irrelevant from the point of view of the selectivity of the model. What matters is that they attract the cations into the filter with their charge, and, at the same time, they exert a hard sphere exclusion and make it difficult for the cations to find space in the filter. Therefore, it proved to be sufficient to model the structural groups of the filter as mobile O 1/2− ions that are confined to the filter region. The profiles shown in Fig.2  The bulk vales, D B i , are taken from experiments; they are 1.334 · 10 −9 , 7.92 · 10 −10 , and 2.032 · 10 −9 m 2 s −1 for Na + , Ca 2+ , and Cl − , respectively. The values inside the cylindrical selectivity filter are fitted to two experimental data points: 250 mM symmetric NaCl at 100 mV for Na + , while added 10 mM luminal CaCl 2 at -100 mV for Ca 2+ . The resulting values (1.27·10 −10 and 1.27·10 −11 for Na + and Ca 2+ , respectively) are fixed and never changed. The Cl − value is irrelevant, because Cl − ions do not carry significant current. These values have been obtained for a moderate system size (H = 54 Å, R = 48 Å). Currents, and, therefore, diffusion coefficients fitted to currents can depend on system size (see later).
Here, we show results for a NaCl-CaCl 2 mixture, where there is 250 mM Na + on both sides of the membrane, while there is 4 µM Ca 2+ on the cytosolic (left) and 50 mM Ca 2+ on the luminal (right) side. The voltage is changed from -150 mV to 150 mV with the ground on the right. The current vs. voltage curves are shown  Figure 3: Currents vs. voltage curves as obtained from experiment [21], the model of Gillespie et al. [21,22] obtained from DFT coupled to the NP equation, and from the NP+LEMC method. In the case of the NP+LEMC method, the currents carried by Na + and Ca 2+ are also shown.
in Fig.3. Total currents are shown in black as obtained from experiments (× symbols), Gillespie's NP+DFT calculations (dashed line), and our NP+LEMC calculations (solid line with filled triangles). Agreement with experiment is very good using both models. The slope of the I − U curve is larger for positive voltages than for negative voltages. More details and understanding can be gained from the current curves for the separate ions, Na + and Ca 2+ (blue and red curves, respectively). At positive voltages, the Ca 2+ current is practically zero, so the total current is carried by Na + ions. At negative voltages, on the other hand, both Na + and Ca 2+ ions contribute to the current.
The explanation of this behavior can be drawn from concentration and electrochemical potential profiles shown in Fig.4. Let us discuss the Na + -profiles first (blue curves). The driving force for the Na + ions is the voltage because Na + concentration is the same on the two sides. The difference between Na + currents at -100 and 100 mV voltages, therefore, is the result of the different competition between Na + and Ca 2+ ions at the two voltages. To understand the difference in this competition, we need to understand the behavior of Ca 2+ ions.
At 100 mV, the driving force for Ca 2+ ions is small because the concentration difference balances the electrical potential difference (see the Ca 2+ electrochemical potential profile in the right panel of Fig.4B). This results in a small Ca 2+ current. Ca 2+ concentration is small in the left bulk, therefore, a Ca 2+ depletion zone is formed in the left vestibule and in the selectivity filter of the channel (beware the logarithmic scale of the concentration axis). That depletion zone results in a decreased concentration of Ca 2+ ions compared to Na + ions inside the channel.
In the case of -100 mV, on the other hand, the Ca 2+ depletion zone in the left vestibule is absent because Ca 2+ ions arrive from the right and "fill up" the left vestibule. There is a more balanced competition between Ca 2+ and Na + ions in the channel and Ca 2+ ions can use the free-energetic advantage that they have over Na + [22,53]. This means that Ca + ions are favored by the crowded selectivity filter, because they provide twice the charge (compared to Na + ions) to balance the charge of O 1/2− ions while occupying about the same space (their diameters are similar: 1.98 and 1.9 Å for Ca 2+ and Na + , respectively). The finite size of the ions plays a crucial role in the selectivity mechanism. This kind of selectivity could not be produced with the PNP theory. Because Ca 2+ concentration is large in the left vestibule of the channel, it drops quickly to the 4 µM value at the left boundary of the solution domain. This introduces a severe system-size dependence into the calculations in this case that is analyzed in Fig.5. Small Ca 2+ concentration on the left hand side corresponds to a large Debye-length in the double layer at left hand side of the membrane (in the left access region). That double layer should fit into the simulation cell (as it does in the case of a larger cell, see black curves in Fig.5).
If the cell is too small (see red curves in Fig.5), there is not enough space for the Ca 2+ concentration to reach the 4 µM limiting value at the left boundary of the cell. The electrochemical potential cannot reach its limiting value either (see the bottom panel). The NP+LEMC calculation, however, provides a solution in this case too, because the layer near the left outer boundary of the cell takes care of the missing access region in an averaged manner. The concentration has a small value in the layer,  while the electrochemical potential profile drops abruptly (large driving force). An appropriate value of the c i ∇µ i product, therefore, is provided by the self-consistent solution so that the continuity equation is satisfied. This solution, however, is approximate. A large portion of the access region with considerable resistance is taken into account in an averaged way by the layer near the system edge. This introduces an error into the value of the Ca 2+ current that is indicated in Fig.5 for both cases. The good behavior of the current-voltage curves compared to experiments and DFT calculations is due to the fit of the Ca 2+ diffusion coefficient, D F Ca 2+ , inside the filter. That value was fitted for negative voltage at the given system size balancing the system-size error.
This result points out the importance of system size in the case of small ionic concentrations, but it also shows the role of the diffusion coefficient profile in NP+LEMC calculations. In confined geometries, the D i (r) profile, although it has a strong relation to the mobility of ions, is primarily an adjustable parameter that is fitted to experiments (as in the present case) or to MD results [6]. The value of D F i takes into account interactions that are absent in the reduced model or accounts for resistances of regions that are absent in the model.
In the case of the RyR channel, for example, the diffusion coefficient profile includes effects of the parts of the ion channel not included in the model: the real RyR channel is much larger than the 46 Å portion modeled here. That region also tunes the total current, but it is not selective. The selectivity filter and its close neigh- borhood modeled here determines ion selectivity and is able to reproduce complex behavior such as anomalous mole fraction experiments discussed previously by Gillespie [21][22][23] and Boda [4].

Rectifying bipolar nanopores
Ion channels are natural nanopores with stable welldefined structures that are very narrow at their selectivity filters (often below 1 nm in radius) to make them suitably selective. The disadvantages, however, are considerable. The structures are often unknown. They are difficult to handle experimentally. Their manipulation is cumbersome with point mutations. Synthetic nanopores, therefore, quickly gained attention due to the fact that they have special properties compared to those of micropores. These special properties arise because the screening length of the electrolyte is comparable to the radius of the nanopore. This fact gives nanopores properties that resemble those of ion channels.
One advantage of synthetic nanopores is that are relatively easy to fabricate [56][57][58][59][60][61][62][63]. They are either solid state nanopores using ion-beam or electron-beam sculpting in, for example, silicon compound membranes, or they are track etched into polymer membranes. Two basic properties of such nanopores are their geometry (shape) and the pattern of surface charge on the pore wall. In our previous works, we studied the bipolar nanopore, where the surface charge is positive on the left hand side of the pore, while it is negative on the right hand side [6,7]. These nanopores rectify ionic current, meaning that they let a much larger amount of ions through at a given sign of voltage than at the opposite sign. In those papers, we used a cylindrical geometry for the nanopore and focused on the effect of the charge pattern, pore radius, and concentration.
Here, we discuss the effect of pore geometry on the rectification properties of a bipolar nanopore. We performed NP+LEMC calculations for three different geometries, the cylindrical (Cyl), single conical (SC), and double conical (DC) shown in Fig.6. Simulations for different voltages have been performed from -150 mV to 150 mV. The concentration of the electrolyte was 0.1 M on both sides. The electrolyte is a 1:1 system (let us call it NaCl), but the diameters of the cations and anions are the same in order to avoid effects from ion size asymmetry. For the same reason, the same diffusion coefficients were used for the two ions.
This makes it possible to focus on the balance of charge and geometrical asymmetries. If the nanopore's shape is symmetric, for example, the cations and anions carry the same amount of current (Cyl and DC). Figure 7 shows the current-voltage relations. Rectification is observed in all three geometries: current is much larger at positive than at negative voltages. The rectifi- cation is defined as the ratio of currents at positive (ON state) and negative voltages (OFF state). The results are shown in Fig.8.
The basic explanation of rectification can be depicted from the concentration profiles (Fig.9). There are depletion zones for cations in the positively charged half region (left), while there are depletion zones for anions in the negatively charged half region (right). In the OFF state (negative voltage) these depletion zones are deeper than in the ON state (positive voltage). More detailed explanation have been given in our previous works [6,7]. Here we focus our discussion to the effect of pore shape.
Total currents are larger in the SC and DC geometries due to their wider entrances. Interestingly, total currents are the same in the SC and DC geometries despite the quite different pore shapes. Whether this is a coincidence or it has a deeper explanation requires further investigation.
What is different in the SC and DC geometries is the partitioning of the total current between Na + and Cl − . While Na + and Cl − currents are the same in the DC geometry due to the symmetric shape, they are different in the SC geometry. The current carried by Na + ions is smaller than the current carried by Cl − ions (middle panel of Fig.7). This is reflected in concentration profiles (see middle panel of Fig.9). The explanation is that the tip of the SC nanopore (left entrance) is positively charged so the depletion zones of Na + ions there is deeper than the depletion zone of Cl − ions on the other side where the pore is wider. Deeper depletion zones result in smaller currents. Therefore, Na + current is smaller in the whole voltage range, but more so in the OFF state at negative voltages.
Rectification is largest in the Cyl geometry because the depletion zones are the deepest in that geometry for both ions. The deepest points of the depletion zones are formed around |z| = 10 Å. There are different effects that form the concentration profiles inside the pore. In the left region, for example, (1) the positive surface charge on the pore wall repels Na + ions, (2) the negative surface charge on the right hand side attracts them, (3) the 0.1 M bulk region acts as a source for the ions, (4) applied potential in the OFF state drives Na + ion to the right, where they have a peak, and (5) Cl − ions (that have a peak on the left) attract them. The balance of all these effects forms the deep depletion zone of Na + at z ≈ −10 Å in the OFF state. Because the Cyl geometry has the smallest radius at the |z| ≈ 10 Å positions, this geometry provides the deepest depletion zone as a result of the dominant effect from the above list, the effect of repelling surface charge.

Summary
We presented the NP+LEMC method that is a hybrid technique, harvesting the advantageous properties of both the NP transport equation (fast calculation of flux) and LEMC particle simulation method (correct calculation of ionic correlations). We applied the method to compute ionic currents through reduced models of an ionic channel and a bipolar nanopore. Reduced models have the advantage of fast calculation (LEMC would not be feasible for an explicit-water model) and intellectual focus. With these models, we can concentrate on those properties of the system that are important to reproduce its behavior as a device.
For the RyR calcium channel, for example, the reduced model is able to reproduce complex selectivity behavior in agreement with experiments by modeling only the "important" amino acids in a reduced way [21,22]. For the bipolar nanopore, the reduced model using implicit water is able to reproduce MD results for an explicit-water model [6]. With further methodological and model development, we intend to simulate nanodevices as close to their real size as possible.