Untitled HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY Vol. 45(1) pp. 73-84 (2017) hjic.mk.uni-pannon.hu DOI: 10.1515/hjic-2017-0011 SIMULATING ION TRANSPORT WITH THE NP+LEMC METHOD. APPLICA- TIONS TO ION CHANNELS AND NANOPORES. DÁVID FERTIG ∗1 , ESZTER MÁDAI1 , MÓNIKA VALISKÓ1 , AND DEZSŐ BODA1,2 1Department of Physical Chemistry, University of Pannonia, Egyetem u. 10., Veszprém, H-8200, HUNGARY 2Institute of Advanced Studies Kőszeg (iASK), Chernel u. 14., Kőszeg, H-9730, HUNGARY 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 es- tablish 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 al- lows 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. Keywords: ion transport, Nernst-Planck, Monte Carlo simulation, nanopore, ion channel 1. Introduction We dedicate this paper to honoring Professor János Liszi, who was the supervisor of one of the authors (DB) and re- spected 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 sim- ulation technique for computing steady-state ion trans- port 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 var- ious problems in the last couple of years. These prob- lems include particle transport through model membranes [1, 2], ion channels [3–5], and nanopores [6–8]. The ba- sic 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 gra- dient 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 advan- tage of including those degrees of freedom of the many- particle 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 chan- nel 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 reduc- tion in the degrees of freedom, however, is that we model water as a dielectric continuum. The effect of using im- plicit solvent instead of an explicit solvent model was dis- cussed 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–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–18] also uses the NP equation to compute flux, but uses a mean-field approach, the Poisson-Boltzmann 74 FERTIG, MÁDAI, VALISKÓ, AND BODA (PB) theory, to describe the statistical mechanical rela- tionship 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–23]. DFT is a con- tinuum 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 the- ory [24]. Its accuracy was demonstrated by computing electrical double layers in extreme conditions and com- paring to Monte Carlo (MC) simulations [25, 26]. The disadvantage of this method is that it can be used effi- ciently only in one dimension. This reduction in dimen- sionality, however, is often feasible provided that we re- duce the model in an intelligent way. 2. Modeling and methodology A good model should be able to explain complex phe- nomena 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: uij(r) =    ∞ if r < Ri + Rj 1 4π�0� qiqj r if r ≥ Ri + Rj (1) where Ri and Rj are the radii, qi and qj are the charges of the different ionic species, �0 is the permittivity of vac- uum, � = 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 rep- resented by two response functions. Energetically, sol- vent acts as a dielectric background that screens the elec- tric field of ions by just dividing by � in the Coulomb- potential (Eq.(1)). Dynamically, solvent molecules col- lide with ions and impede their diffusion; we take that effect into account by a diffusion coefficient function, Di(r). The ions of the electrolyte diffuse through a pore be- tween two bulk containers (see Fig.1). The pore pen- etrates a membrane that separates the two containers. The pore and the membrane are also modeled in a re- duced way; they are confined by hard walls with which the hard sphere ion cannot overlap. Other details of the Figure 1: Simulation cell used in the application of the NP+LEMC method. The model is rotationally symmetric: 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 poten- tial in this region is not constant, but changes from one constant value to another (the values in the two bulk re- gions) in a monotonic manner to provide the driving force of ionic transport. Parameters R and H characterize sys- tem size. 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; concen- trations 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 drift- diffusion. Our procedure is a hybrid method that sepa- rates the configuration degrees of freedom of ions (ion positions) from the kinetic degrees of freedom (ion ve- locities). 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 cur- rent 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 peculiari- ties beyond ideal-gas behavior are from the potential en- ergy, that, in turn, depends only on the configuration coor- dinates [27–30]. Out of equilibrium, however, this sepa- ration 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 con- Hungarian Journal of Industry and Chemistry SIMULATING ION TRANSPORT WITH NP+LEMC 75 tinuum 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 equilib- rium. There is no well-established non-equilibrium sta- tistical mechanics finding its way into textbooks despite the fact that many authors made considerable advances in this field [31–36]. The most problematic quantities are those containing entropic effects, such as the chem- ical 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 ho- mogeneity defines thermodynamic equilibrium (let us as- sume 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 (electrod- iffusion) is the NP equation, ji(r) = − 1 kT Di(r)ci(r)∇µi(r), (2) where Di(r) is the diffusion coefficient profile of ionic species i, ci(r) is the concentration profile, µi(r) is the electrochemical potential profile, ji(r) is the particle flux density, k is Boltzmann’s constant and T = 298.15K is the temperature. The resulting ji(r) must satisfy the continuity equation: ∇ · ji(r) = 0. (3) 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 statisti- cal mechanical task to averaging over states in the con- figuration space just as we did in the case of equilibrium statistical mechanics. If we treat Di(r) as a parameter, the task is reduced to finding the proper relation between ci(r) and µi(r). Defining local concentration in a non- equilibrium situation is the same as in equilibrium: we compute the average number of particles in a small vol- ume 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 lo- cal equilibrium (LE). If we divide the simulation cell into subvolumes Bα, we can characterize this subvol- ume 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 el- ements. The assumption of LE is also present in other methods (though not necessarily stated explicitly), where the ci(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 sub- volumes are open systems and they are in LE with a con- stant volume (V α), temperature (T ), and electrochemical potential (µαi ). Thus, we suggested to use Grand Canoni- cal Monte Carlo (GCMC) simulations, where particle in- sertion/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): min(1; pαi,χ(r)), where pαi,χ(r) = Nαi !(V α)χ (Nαi + χ)! exp ( − ∆U(r) − χµαi kT ) . (4) 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 po- sition 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 ac- count 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 config- urations 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 el- ements and particle species, respectively. For compari- son, 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 conserva- tion of mass (∇ · ji(r) = 0) is satisfied. The procedure can be summarized as µαi [n] LEMC −−−−→ cαi [n] NP −−→ jαi [n] ∇·j=0 −−−−→ µαi [n + 1]. (5) The electrochemical potentials for the next iteration, µαi [n + 1], are computed from the results of the previous iteration, cαi [n], on the basis of the divergence-theorem (also known as Gauss-Ostrogradsky’s theorem). The con- tinuity equation is converted to a surface integral: 0 = ∫ Bα ∇ · ji(r) dV = ∮ Sα ji(r) · n(r) da, (6) where volume Bα is bounded by surface Sα and n(r) denotes the normal vector pointing outward at position 45(1) pp. 73-84 (2017) 76 FERTIG, MÁDAI, VALISKÓ, AND BODA 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 , and D̂ αβ 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 consti- tute Sα: 0 = ∑ β,Sαβ∈Sα ĵ αβ i · n αβaαβ, (7) 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. An appropriately chosen initial set of electrochem- ical potentials is chosen (µαi [1]; in general: µ α i [n], where [n] denotes the nth iteration). 2. Using these µαi [n] parameters as inputs, LEMC sim- ulations are performed. The resulting concentrations are denoted by cαi [n]. 3. The flux computed from the {µαi [n], c α i [n]} pair usu- ally does not satisfy Eq.(7). The next set of electro- chemical potential is calculated by assuming that the {µαi [n + 1], c α i [n]} pair does satisfy Eq.(7). If we write the value of ĵαβi as given by the NP equation into Eq.(7), we obtain 0 = ∑ β D̂ αβ i ĉ αβ i [n]∇µ̂ αβ i [n + 1] · n αβ aαβ, (8) where β runs over all the surface elements Sαβ that constitute Sα. Eq.(8) is a system of linear equations; the unknown variables are denoted by µα,CALi [n + 1], where CAL refers to the fact that these values come from calculations by solving Eq.(8). 4. To achive faster and more robust convergence in the case of large driving forces, the electrochemical po- tential 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: µ α,MIX i [n + 1] = biµ α,CAL i [n + 1] + (1 − bi)µ α,MIX i [n], (9) where bi 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 bi values prevent this fluctuation at the price of making the convergence slower. 5. The input of the (n + 1)th LEMC simulation is µ α,MIX i [n + 1]. The system of linear equations contains the boundary conditions for the electrochemical potential and concen- tration 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 µLi = kT ln c L i + µ EX,L i + qiΦ L (10) and µRi = kT ln c R i + µ EX,R i + qiΦ R, (11) where µEX,Li and µ EX,R i are the excess chemical poten- tials in the absence of an external field (determined by the Adaptive GCMC method of Malasics et al. [38]), while qiΦ L and qiΦ R are the interactions with the applied elec- trical 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 concen- trations 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 interac- tions between particles (and interactions of particles with the pore), but also the interaction with an external elec- trical potential, Φappl(r). This applied potential is calcu- lated by solving Laplace’s equation ∇2Φappl(r) = 0 (12) 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 z- axis 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 disad- vantageous 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 po- tential (containing also the interaction with the applied field) is µαi − kT ln c α i and can be computed as exp [− (µαi /kT − ln c α i )] = 〈exp [−∆U α i /kT ]〉 , (13) where ∆Uαi is the energy change associated with the in- sertion of a test particle of species i into a randomly cho- sen position in the elementary cell α. This is the same en- ergy computed in the particle insertion step of the LEMC technique, therefore, it does not require additional com- putational cost. The concentration can be calculated from Hungarian Journal of Industry and Chemistry SIMULATING ION TRANSPORT WITH NP+LEMC 77 Eq.(13) because µαi is known (that is the input of the sim- ulation) 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 sam- pling 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 bet- ter, when concentrations are very large in the volume el- ement. 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 solu- tion 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 av- eraging: Ji(z) = 2π ∫ R(z) 0 rji(z, r) · nz dr (14) where R(z) denotes the radius of the pore at coordinate z (|z| < Hmemb/2, where Hmemb is the thickness of the membrane) and nz is the unit vector along the z-axis. The electrical current of an ion is then Ii = qiJi, (15) and the total current is I = ∑ i qiJi. (16) 3. Results and discussion 3.1 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 Ca2+ selective channels, meaning that the presence of only a µM Ca2+ in the bulk decreases the Na+ current to half of its value compared to its value in the absence of Ca2+. This is called micromolar Ca2+ affinity (or Ca2+ block) because a very small amount of Ca2+ is enough to block the channel [44]. The price of high Ca2+ selectivity is low Ca2+ current through this channel. The Ca2+ ions flow into the muscle cell when the L-type calcium channels open in response to an electrical signal (action potential) in a small quan- tity that is not sufficient to initiate muscle contraction. The solution of Nature for this problem is an amplifica- tion mechanism. A large amount of Ca2+ 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 mem- brane of the SR that are opened by the Ca2+ ions pro- vided 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 Ca2+ flux that is necessary for muscle contraction by binding to tropomyosins on the actin filaments and al- lowing myosin to climb up the filament. These channels are wider (also, less Ca2+ selective) than the L-type channels. A large amount of experimen- tal data is available for this channel [45–47]. On the basis of this database, Gillespie et al. developed a model for the RyR channel [21, 22]. Although the model was a one- dimensional reduction, studied with a one-dimensional NP+DFT, they were able to reproduce hundreds of ex- perimental current-voltage curves. Later, when the NP+LEMC method became avail- able, 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–54] for the L-type cal- cium 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 im- portant part of the channel is the selectivity filter, the bot- tleneck 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 lo- cus). 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–55] that the strong Ca2+ selectivity can be reproduced by a reduced model where eight half-charged oxygen atoms (O1/2−) of the four COO− groups force a strict competition between Ca2+ and Na+ ions for space in the crowded filter. The exact positions of the O1/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 ex- ert 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 O1/2− ions that are confined to the filter re- gion. The profiles shown in Fig.2 are results of the sim- ulations and show how the eight O1/2− ions distribute in a given segment of the channel. There are 32 such O1/2− ions, plus point charges placed on a ring on the luminal side. The geometrical features (filter radius is 5.5 Å, filter length is 15 Å, and radii of left and right vestibules are 22 and 9 Å, respectively) are also designed on the basis of Gillespie’s model that, in turn, are based on fitting to experimental data. This model was designed and used unchanged in ev- 45(1) pp. 73-84 (2017) 78 FERTIG, MÁDAI, VALISKÓ, AND BODA Figure 2: The three-dimensional model [4] of the RyR calcium channel based on the model of Gillespie et al. [21, 22]. Each curve shows the distribution of eight O1/2− ions that are confined to the given region but otherwise free to move there. The charges of the E4902 amino acids are represented as point charges (−1/2e) at fixed positions on a ring. The thickness of the membrane is Hmemb = 46 Å. The dielectric constant is � = 78.5. ery calculation. What remains to be specified are the dif- fusion coefficient profiles of the various ionic species. These profiles depend only on z, Di(z); we assume no variation in the radial dimension. The values outside the channel and in the selectivity filter, DBi and D F i , respec- tively, are constant. In the vestibules at the two entrances of the channel, the diffusion constant profiles are interpo- lated between the DBi and D F i values in a way that their value changes in proportion with the cross section of the channel. The bulk vales, DBi , are taken from experiments; they are 1.334·10−9, 7.92·10−10, and 2.032·10−9 m2s−1 for Na+, Ca2+, 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 CaCl2 at -100 mV for Ca 2+. The resulting values (1.27·10−10 and 1.27·10−11 for Na+ and Ca2+, 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-CaCl2 mixture, where there is 250 mM Na+ on both sides of the mem- brane, while there is 4 µM Ca2+ on the cytosolic (left) and 50 mM Ca2+ on the luminal (right) side. The volt- age 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 ex- periment [21], the model of Gillespie et al. [21, 22] ob- tained 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 Ca2+ are also shown. in Fig.3. Total currents are shown in black as obtained from experiments (× symbols), Gillespie’s NP+DFT cal- culations (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 posi- tive voltages than for negative voltages. More details and understanding can be gained from the current curves for the separate ions, Na+ and Ca2+ (blue and red curves, respectively). At positive voltages, the Ca2+ current is practically zero, so the total current is carried by Na+ ions. At negative voltages, on the other hand, both Na+ and Ca2+ 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 competi- tion between Na+ and Ca2+ ions at the two voltages. To understand the difference in this competition, we need to understand the behavior of Ca2+ ions. At 100 mV, the driving force for Ca2+ ions is small because the concentration difference balances the electri- cal potential difference (see the Ca2+ electrochemical po- tential profile in the right panel of Fig.4B). This results in a small Ca2+ current. Ca2+ concentration is small in the left bulk, therefore, a Ca2+ 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 Ca2+ ions compared to Na+ ions inside the channel. In the case of -100 mV, on the other hand, the Ca2+ depletion zone in the left vestibule is absent because Ca2+ ions arrive from the right and “fill up” the left Hungarian Journal of Industry and Chemistry SIMULATING ION TRANSPORT WITH NP+LEMC 79 (A) (B) Figure 4: (A) Concentration and (B) electrochemical po- tential profiles of Na+ and Ca2+ for two opposing volt- ages: -100 mV (left panels) and 100 mV (right pan- els). The striped parts indicate the cytosolic and luminal vestibules, while the gray area is the selectivity filter. vestibule. There is a more balanced competition between Ca2+ and Na+ ions in the channel and Ca2+ 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 O1/2− ions while occupying about the same space (their diameters are similar: 1.98 and 1.9 Å for Ca2+ 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 Ca2+ 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 in- troduces a severe system-size dependence into the calcu- lations in this case that is analyzed in Fig.5. Small Ca2+ 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 Ca2+ 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 calcula- tion, however, provides a solution in this case too, be- cause 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, Figure 5: Concentration (top panel) and electrochemical potential (bottom panel) profiles of Ca2+ for -100 mV. The results of simulations for two system sizes are shown: H = 54 Å (red) and H = 180 Å (black). while the electrochemical potential profile drops abruptly (large driving force). An appropriate value of the ci∇µi product, therefore, is provided by the self-consistent so- lution so that the continuity equation is satisfied. This solution, however, is approximate. A large por- tion 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 Ca2+ 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 Ca2+ diffusion coefficient, DF Ca2+ , 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 Di(r) profile, although it has a strong relation to the mobility of ions, is primarily an adjustable parameter that is fitted to experi- ments (as in the present case) or to MD results [6]. The value of DFi takes into account interactions that are ab- sent 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 dif- fusion 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- 45(1) pp. 73-84 (2017) 80 FERTIG, MÁDAI, VALISKÓ, AND BODA Figure 6: Bipolar nanopore geometries: Cyl: cylindrical; SC: single conical; DC: double conical. The pore is pos- itively charged on the left hand side (−30 < z < 0 Å), while negatively charged on the right hand side (0 < z < 30 Å) keeping the surface charge fixed (σ = ±0.1 C/m2). Minimal and maximal pore radii are 10 and 20 Å, respec- tively. borhood modeled here determines ion selectivity and is able to reproduce complex behavior such as anomalous mole fraction experiments discussed previously by Gille- spie [21–23] and Boda [4]. 3.2 Rectifying bipolar nanopores Ion channels are natural nanopores with stable well- defined structures that are very narrow at their selectivity filters (often below 1 nm in radius) to make them suit- ably selective. The disadvantages, however, are consider- able. The structures are often unknown. They are difficult to handle experimentally. Their manipulation is cumber- some with point mutations. Synthetic nanopores, therefore, quickly gained atten- tion due to the fact that they have special properties com- pared 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 rela- tively easy to fabricate [56–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 proper- ties of such nanopores are their geometry (shape) and the pattern of surface charge on the pore wall. Figure 7: Current-voltage relations for the three nanopore geometries (Cyl, SC, and DC from bottom to top). Cur- rents carried by Na+ (solid blue), Cl− (dashed red), and their sum (dot-dashed black) are shown as a function of voltage. The insets magnify the results for negative volt- ages. 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, mean- ing 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 ra- dius, and concentration. Here, we discuss the effect of pore geometry on the rectification properties of a bipolar nanopore. We per- formed NP+LEMC calculations for three different ge- ometries, the cylindrical (Cyl), single conical (SC), and double conical (DC) shown in Fig.6. Simulations for dif- ferent 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. Rectifi- cation is observed in all three geometries: current is much larger at positive than at negative voltages. The rectifi- Hungarian Journal of Industry and Chemistry SIMULATING ION TRANSPORT WITH NP+LEMC 81 Figure 8: Rectification (defined as the absolute value of the current ratio in the ON and OFF states) as a function of the absolute value of the voltage for the three nanopore geometries. 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 deple- tion 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 expla- nation 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 investiga- tion. 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 differ- ent 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 Figure 9: Concentration profiles of Na+ (blue) and Cl− (red) in the ON (solid lines) and OFF (dashed lines) states for the three nanopore geometries (Cyl, SC, and DC from bottom to top). 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. 4. 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 chan- nel and a bipolar nanopore. Reduced models have the ad- vantage 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 re- duced model is able to reproduce complex selectivity behavior in agreement with experiments by modeling only the “important” amino acids in a reduced way 45(1) pp. 73-84 (2017) 82 FERTIG, MÁDAI, VALISKÓ, AND BODA [21, 22]. For the bipolar nanopore, the reduced model us- ing implicit water is able to reproduce MD results for an explicit-water model [6]. With further methodological and model development, we intend to simulate nanode- vices as close to their real size as possible. Acknowledgement The financial support of the National Research, Develop- ment and Innovation Office (NKFIH K124353) is grate- fully acknowledged. The present article was published in the frame of the Project No. GINOP-2.3.2-15-2016- 00053 and supported by the UNKP-17-4 (to M.V.) and UNKP-17-1 (to E.M.), the New National Excellence Pro- gram of the Ministry of Human Capacities. We thank Dirk Gillespie for helpful discussions. REFERENCES [1] Boda, D., Gillespie, D.: Steady state electrodif- fusion from the Nernst-Planck equation coupled to Local Equilibrium Monte Carlo simulations, J. Chem. Theor. Comput., 2012 8(3), 824–829, DOI: 10.1021/ct2007988 [2] Ható, Z., Boda, D., Kristóf, T.: Simulation of steady-state diffusion: Driving force ensured by Dual Control Volumes or Local Equilibrium Monte Carlo, J. Chem. Phys., 2012 137(5), 054109, DOI: 10.1063/1.4739255 [3] Boda, D., Kovács, R., Gillespie, D., Kristóf, T.: Se- lective transport through a model calcium channel studied by Local Equilibrium Monte Carlo simula- tions coupled to the Nernst-Planck equation, J. Mol. Liq., 2014 189, 100, DOI: 10.1016/j.molliq.2013.03.015 [4] Boda, D.: Monte Carlo Simulation of Electrolyte Solutions in Biology: In and Out of Equilib- rium, Annual Reports in Computational Chemistry, 2014, vol. 10, chap. 5 (Elsevier), 127–163, DOI: 10.1016/b978-0-444-63378-1.00005-7 [5] Ható, Z., Boda, D., Gillepie, D., Vrabec, J., Rutkai, G., Kristóf, T.: Simulation study of a rectifying bipolar ion channel: detailed model versus reduced model, Cond. Matt. Phys., 2016 19(1), 13802, DOI: 10.5488/cmp.19.13802 [6] Ható, Z., Valiskó, M., Kristóf, T., Gillespie, D., Boda, D.: Multiscale modeling of a rectifying bipolar nanopore: explicit-water versus implicit- water simulations, Phys. Chem. Chem. Phys., 2017 17(27), 17816–17826, DOI: 10.1039/C7CP01819C [7] Matejczyk, B., Valiskó, M., Wolfram, M.T., Pietschmann, J.F., Boda, D.: Multiscale modeling of a rectifying bipolar nanopore: Comparing Poisson- Nernst-Planck to Monte Carlo, J. Chem. Phys., 2017 146(12), 124125, DOI: 10.1063/1.4978942 [8] Mádai, E., Valiskó, M., Dallos, A., Boda, D.: Sim- ulation of a model nanopore sensor: Ion competi- tion underlies device behavior, J. Chem. Phys., 2017 147(24), 244702, DOI: 10.1063/1.5007654 [9] van Gunsteren, W., Berendsen, H.: Algorithms for Brownian dynamics, Mol. Phys., 1982 45(3), 637– 647, DOI: 10.1080/00268978200100491 [10] Turq, P., Lantelme, F., Friedman, H.L.: Brown- ian dynamics: Its application to ionic solutions, J. Chem. Phys., 1977 66(7), 3039–3044, DOI: 10.1063/1.434317 [11] Corry, B., Kuyucak, S., Chung, S.H.: Tests of con- tinuum theories as models of ion channels. II. Poisson-Nernst-Planck theory versus Brownian dy- namics, Biophys. J., 2000 78(5), 2364–2381, DOI: 10.1016/s0006-3495(00)76781-6 [12] Chung, S.H., Allen, T.W., Hoyles, M., Kuyucak, S.: Permeation Of Ions Across The Potassium Chan- nel: Brownian Dynamics Studies, Biophys. J., 1999 77(5), 2517–2533, DOI: 10.1016/s0006-3495/99/77087-6 [13] Chung, S.H., Hoyles, M., Allen, T., Kuyucak, S.: Study Of Ionic Currents Across A Model Mem- brane Channel Using Brownian Dynamics, Bio- phys. J., 1998 75(2), 793–809, DOI: 10.1016/s0006- 3495(98)77569-1 [14] Berti, C., Furini, S., Gillespie, D., Boda, D., Eisen- berg, R.S., Sangiorgi, E., Fiegna, C.: A 3-D Brow- nian Dynamics simulator for the study of ion permeation through membrane pores, J. Chem. Theor. Comput., 2014 10(8), 2911–2926, DOI: 10.1021/ct4011008 [15] Zheng, Q., Chen, D., Wei, G.W.: Second-order Poisson-Nernst-Planck solver for ion transport, J. Comp. Phys., 2011 230(13), 5239–5262, DOI: 10.1016/j.jcp.2011.03.020 [16] Schuss, Z., Nadler, B., Eisenberg, R.S.: Deriva- tion of Poisson and Nernst-Planck equations in a bath and channel from a molecular model, Phys. Rev. E, 2001 6403(3), 036116, DOI: 10.1103/phys- reve.64.036116 [17] Nonner, W., Eisenberg, B.: Ion permeation and glu- tamate residues linked by Poisson-Nernst-Planck theory in L-type calcium channels, Biophys. J., 1998 75(3), 1287–1305, DOI: 10.1016/s0006- 3495(98)74048-2 [18] Pietschmann, J.F., Wolfram, M.T., Burger, M., Trautmann, C., Nguyen, G., Pevarnik, M., Bayer, V., Siwy, Z.: Rectification properties of conically shaped nanopores: consequences of miniaturiza- tion, Phys. Chem. Chem. Phys., 2013 15(3), 16917– 16926, DOI: 10.1039/C3CP53105H [19] Gillespie, D., Nonner, W., Eisenberg, R.S.: Cou- pling Poisson-Nernst-Planck and density functional theory to calculate ion flux, J. Phys.: Cond. Matt., 2002 14(46), 12129–12145, DOI: 10.1088/0953- 8984/14/46/317 [20] Gillespie, D., Nonner, W., Eisenberg, R.S.: Density functional theory of charged, hard-sphere fluids, Phys. Rev. E, 2003 68(3), 031503, DOI: 10.1103/phys- reve.68.031503 [21] Gillespie, D., Xu, L., Wang, Y., Meissner, G.: (De)constructing the ryanodine receptor: Modeling Hungarian Journal of Industry and Chemistry SIMULATING ION TRANSPORT WITH NP+LEMC 83 ion permeation and selectivity of the calcium re- lease channel, J. Phys. Chem. B, 2005 109(32), 15598–15610, DOI: 10.1021/jp052471j [22] Gillespie, D.: Energetics of Divalent Selectivity in a Calcium Channel: The Ryanodine Receptor Case Study, Biophys. J., 2008 94(4), 1169–1184, DOI: 10.1529/biophysj.107.116798 [23] Gillespie, D., Giri, J., Fill, M.: Reinterpreting the anomalous mole fraction effect: The Ryanodine re- ceptor case study, Biophys. J., 2009 97(8), 2212– 2221, DOI: 10.1016/j.bpj.2009.08.009 [24] Rosenfeld, Y.: Free-energy model for the inho- mogeneous hard-sphere fluid mixture and density- functional theory of freezing, Phys. Rev. Lett., 1989 63(9), 980–983, DOI: 10.1103/physrevlett.63.980 [25] Gillespie, D., Valiskó, M., Boda, D.: Density func- tional theory of the electrical double layer: the RFD functional, J. Phys.-Cond. Matt., 2005 17(42), 6609–6626, DOI: 10.1088/0953-8984/17/42/002 [26] Valiskó, M., Gillespie, D., Boda, D.: Selective ad- sorption of ions with different diameter and valence at highly-charged interfaces, J. Phys. Chem. C, 2007 111(43), 15575–15585, DOI: 10.1021/jp073703c [27] McQuarrie, D.A.: Statistical Mechanics (Univer- sity Science Books, Sausalito), 2000, ISBN: 9781891389153 [28] Allen, M.P., Tildesley, D.J.: Computer Simulation of Liquids, Oxford Science Publ (Clarendon Press, New York), 1989, ISBN: 9780198556459 [29] Frenkel, D., Smit, B.: Understanding molecular simulations (Academic Press, San Diego), 1996, ISBN: 9780122673702 [30] Hansen, J.P., McDonald, I.R.: Theory of Simple Liquids (Elsevier, Academic Press), 2006 ISBN: 9780080455075 [31] Prigogine, I.: Introduction to Thermodynamics of Irreversible Processes (John Wiley and Sons, New York), 1968, ISBN: 9780470699287 [32] Kreuzer, H.J.: Nonequilibrium Thermodynamics and its Statistical Foundations, Monographs on the physics and chemistry of materials (Clarendon Press, Oxford), 1983, ISBN: 9780198513759 [33] de Groot, S.R., Mazur, P.: Non-Equilibrium Ther- modynamics, Dover Books on Physics (Dover, New York), 2013, ISBN: 9780486153506 [34] Tschoegl, N.W.: Fundamentals of Equilibrium and Steady-State Thermodynamics (Elsevier, Amster- dam), 2000, ISBN: 9780080532110 [35] Attard, P.: Thermodynamics and Statistical Mechanics, Equilibrium by Entropy Maximi- sation (Academic Press, London), 2002, ISBN: 9780120663217 [36] Evans, D.J., Morriss, G.: Statistical Mechanics of Nonequilibrium Liquids (Cambridge University Press, New York), 2008, ISBN: 9781139471930 [37] Beck, T.L., Paulaitis, M.E., Pratt, L.R.: The Poten- tial Distribution Theorem and Models of Molecu- lar Solutions (Cambridge University Press, Cam- bridge), 2006, ISBN: 9781139457637 [38] Malasics, A., Boda, D.: An efficient iterative grand canonical Monte Carlo algorithm to determine in- dividual ionic chemical potentials in electrolytes, J. Chem. Phys., 2010 132(24), 244103, DOI: 10.1063/1.3443558 [39] Boda, D., Gillespie, D., Nonner, W., Henderson, D., Eisenberg, B.: Computing induced charges in inhomogeneous dielectric media: Application in a Monte Carlo simulation of complex ionic systems, Phys. Rev. E, 2004 69(4), 046702, DOI: 10.1103/phys- reve.69.046702 [40] Boda, D., Valiskó, M., Eisenberg, B., Nonner, W., Henderson, D., Gillespie, D.: The effect of protein dielectric coefficient on the ionic selectivity of a cal- cium channel, J. Chem. Phys., 2006 125(3), 034901, DOI: 10.1063/1.2212423 [41] Widom, B.: Some Topics in the Theory of Flu- ids, J. Chem. Phys., 1963 39(11), 2808–2812, DOI: 10.1063/1.1734110 [42] Widom, B.: Structure of interfaces from uniformity of the chemical potential, J. Stat. Phys., 1978 19(6), 563–574, DOI: 10.1007/bf01011768 [43] Hille, B.: Ion channels of excitable membranes (Sinauer Associates, Sunderland), 2001, ISBN: 9780878933211 [44] Sather, W.A., McCleskey, E.W.: Permeation and selectivity in calcium channels, Ann. Rev. Phys- iology, 2003 65(1), 133–159, DOI: 10.1146/an- nurev.physiol.65.092101.142345 [45] Gao, L., Balshaw, D., Xu, L., Tripathy, A., Xin, C., Meissner, G.: Evidence for a Role of the Lume- nal M3-M4 Loop in Skeletal Muscle Ca2+ Release Channel (Ryanodine Receptor) Activity and Con- ductance, 2000 79(2), 828–840, DOI: 10.1016/s0006- 3495(00)76339-9 [46] Wang, Y., Xu, L., Pasek, D.A., Gillespie, D., Meiss- ner, G.: Probing the Role of Negatively Charged Amino Acid Residues in Ion Permeation of Skeletal Muscle Ryanodine Receptor, 2005 89(1), 256–265, DOI: 10.1529/biophysj.104.056002 [47] Xu, L., Wang, Y., Gillespie, D., Meissner, G.: Two Rings of Negative Charges in the Cytosolic Vestibule of Type-1 Ryanodine Receptor Modulate Ion Fluxes, 2006 90(2), 443–453, DOI: 10.1529/bio- physj.105.072538 [48] Boda, D., Valiskó, M., Eisenberg, B., Nonner, W., Henderson, D., Gillespie, D.: Combined ef- fect of pore radius and protein dielectric coeffi- cient on the selectivity of a calcium channel, Phys. Rev. Lett., 2007 98(16), 168102, DOI: 10.1103/phys- revlett.98.168102 [49] Gillespie, D., Boda, D.: The anomalous mole frac- tion effect in calcium channels: A measure of pref- erential selectivity, Biophys. J., 2008 95(6), 2658– 2672, DOI: 10.1529/biophysj.107.127977 45(1) pp. 73-84 (2017) 84 FERTIG, MÁDAI, VALISKÓ, AND BODA [50] Boda, D., Valiskó, M., Henderson, D., Eisenberg, B., Gillespie, D., Nonner, W.: Ion selectivity in L- type calcium channels by electrostatics and hard- core repulsion, J. Gen. Physiol., 2009 133(5), 497– 509, DOI: 10.1085/jgp.200910211 [51] Malasics, M., Boda, D., Valiskó, M., Henderson, D., Gillespie, D.: Simulations of calcium channel block by trivalent ions: Gd3+ competes with per- meant ions for the selectivity filter, Biochim. et Bio- phys. Acta - Biomembranes, 2010 1798(11), 2013– 2021, DOI: 10.1016/j.bbamem.2010.08.001 [52] Rutkai, G., Boda, D., Kristóf, T.: Relating bind- ing affinity to dynamical selectivity from dynamic Monte Carlo simulations of a model calcium chan- nel, J. Phys. Chem. Lett., 2010 1(14), 2179–2184, DOI: 10.1021/jz100718n [53] Boda, D., Giri, J., Henderson, D., Eisenberg, B., Gillespie, D.: Analyzing the components of the free energy landscape in a calcium selective ion chan- nel by Widom’s particle insertion method, J. Chem. Phys., 2011 134(5), 055102, DOI: 10.1063/1.3532937 [54] Boda, D., Henderson, D., Gillespie, D.: The role of solvation in the binding selectivity of the L-type cal- cium channel, J. Chem. Phys., 2013 139(5), 055103, DOI: 10.1063/1.4817205 [55] Nonner, W., Catacuzzeno, L., Eisenberg, B.: Bind- ing and selectivity in L-type calcium channels: A mean spherical approximation, Biophys. J., 2000 79(4), 1976–1992, DOI: 10.1016/s0006-3495(00)76446-0 [56] Siwy, Z., Fulinski, A.: Fabrication of a synthetic nanopore ion pump, Phys. Rev. Lett., 2002 89(19), 198103, DOI: 10.1103/physrevlett.89.198103 [57] Siwy, Z., Apel, P., Baur, D., Dobrev, D.D., Korchev, Y.E., Neumann, R., Spohr, R., Trautmann, C., Voss, K.O.: Preparation of synthetic nanopores with trans- port properties analogous to biological channels, Surf. Sci., 2003 532, 1061–1066, DOI: 10.1016/s0039- 6028(03)00448-5 [58] Siwy, Z., Apel, P., Dobrev, D., Neumann, R., Spohr, R., Trautmann, C., Voss, K.: Ion transport through asymmetric nanopores prepared by ion track etch- ing, Nuclear Instruments & Methods In Phys. Re- search Section B-beam Interactions With Materi- als Atoms, 2003 208, 143–148, DOI: 10.1016/s0168- 583x(03)00884-x [59] Howorka, S., Siwy, Z.: Nanopores: Generation, En- gineering, and Singly-Molecule Applications, in Handbook of Single-Molecule Biophysics (P. Hin- terdorfer, A. van Oijen, eds.), Advances in Chemi- cal Physics, chap. Chapter 11 (Springer), 2009 293– 339, DOI: 10.1007/978-0-387-76497-9_11 [60] Guo, W., Tian, Y., Jiang, L.: Asymmetric Ion Transport through Ion-Channel-Mimetic Solid- State Nanopores, Acc. Chem. Res., 2013 46(12), 2834–2846, DOI: 10.1021/ar400024p [61] Gibb, T., Ayub, M.: Solid-State Nanopore Fabri- cation, in Engineered Nanopores for Bioanalytical Applications (Elsevier BV), 2013 121–140, DOI: 10.1016/b978-1-4377-3473-7.00005-4 [62] Guan, W., Li, S.X., Reed, M.A.: Voltage gated ion and molecule transport in engineered nanochan- nels: theory, fabrication and applications, Nanotech- nology, 2014 25(12), 122001, DOI: 10.1088/0957- 4484/25/12/122001 [63] Zhang, H., Tian, Y., Jiang, L.: Fundamen- tal studies and practical applications of bio- inspired smart solid-state nanopores and nanochan- nels, Nano Today, 2016 8(3), 1470–1478, DOI: 10.1016/j.nantod.2015.11.001 Hungarian Journal of Industry and Chemistry