Our devices are getting smaller and smaller. Nanodevices that involve ion transport generally have regions that can be considered macroscopic from the molecular point of view, with characteristic size of m scale and larger (e.g., wide part of a conical nanopore, access regions, bulk regions with electrodes). Continuum theories are often sufficient to model and compute devices built of such regions. These reduced (or coarse-grained) models make approximations like point-charge ions and implicit water as a trade-off for fast computation time. As the dimensions of the devices shrink, however, nanoscopic regions appear (e.g., narrow bottlenecks of nanopores), where the assumptions of these reduced models are no longer valid and molecular models and particle simulations like all-atom molecular dynamics (MD) are needed to reveal molecular mechanisms. To connect with experimental data then requires a multiscale modeling approach, where different models/methods with different resolutions are applied to different regions. In our approach, multiscaling is applied between models for the whole system in separate simulations rather than between separate regions in a single model and simulation. The multiscaling part of our procedure is connecting the various modeling levels; molecular insights obtained from atomic models, for example, can be used to design better reduced models.

Unexpectedly, however, reduced models often reproduce—and predict—experimental data, even though they should, in principle, not work because their underlying atomic physics is too approximate; one should not assume that water molecules are absent or that ions are point charges when the pore they go through is only a few times wider than the “real” ion or water. Yet, the Poisson-Nernst-Planck (PNP) theory can reproduce experimental data for rectifying nanopores [1, 2, 3, 4, 5, 6]. Other examples include reduced models using hard sphere ions reproducing and predicting nonlinear phenomena in biological ion channels [7, 8, 9, 10, 11, 12, 13, 14] and in nanopores [15, 16]. Moreover, in many cases including nanopores and ion channels, reduced models are the only way to connect with experimental results, as the required low ion concentrations and small applied voltages are inaccessible to all-atom MD simulations.

The goal of this work is to report computational results for a nanodevice, a rectifying nanopore, on two distinct modeling levels that differ in the treatment of water. In one model, water is a continuum background and the ions are charged, hard spheres that move via drift-diffusion. In the other model, MD simulations with explicit water molecules are used. This simple case study allows us to identify the effects of using an implicit water model versus having corpuscular water molecules.

Consistent with the work cited above, we find that both models produce qualitatively the same input-output relations (e.g., current vs. voltage, pore surface charge, or ion concentration). On the other hand, the molecular-scale picture inside the nanopore for the two models is quite different, as one would expect. The solution to this paradox is that both models produce qualitatively the same axial ion concentration profiles (i.e., cross-section averaged profiles), and because these are the first-order determinants of current, the two models produce the same device characteristics.

Our results suggest that reduced models that reproduce experimental results are, in fact, valid because they capture the overall device physics correctly despite the incomplete description of the molecular level phenomena. Therefore, they are useful for understanding the device-level physics that produces the input-output relations and for device design. At that level, the molecular details are of second order importance. However, to understand the physics at an atomic scale (e.g., ion and water structure inside a confining nanopore), more detailed molecular models must be used.


In an all-atom model, water molecules are modeled explicitly and the system is simulated with the MD technique. The ionic current, driven by an external electric field, is simulated directly by counting the ions that passed through the pore (more detail is found in the Appendix and in the ESI).

Coarse-grained models reduce the number of degrees of freedom in the Hamiltonian by integrating out some of these degrees of freedom and replacing them with a response function, or simplify model potentials describing certain parts of the system. In this work, the reduced model describes water as an implicit continuum background, where the effect of water is replaced by various response functions:
(1) The ability of water to affect the movement of ions through friction (a dynamic effect) is taken into account by a diffusion coefficient, , in the Nernst-Planck (NP) transport equation with which we compute the flux:


where is the particle flux density of ionic species , is concentration, is Boltzmann’s constant, is temperature, and is the electrochemical potential profile. Its gradient is the driving force of the steady-state transport.
(2) The ability of water to screen the charges of ions (an energetic effect) is taken into account by a dielectric constant, , in the denominator of the Coulomb-potential acting between the charged hard spheres with which we model the ions:


where is the charge and is the radius of ionic species , is the permittivity of vacuum, and is the distance between the ions.

We simulate this model with the Local Equilibrium Monte Carlo (LEMC) technique [17, 18, 19], which is practically a grand canonical Monte Carlo simulation devised for a non-equilibrium situation. The input variable of the LEMC simulation is the chemical potential profile, , that is not constant for a system out of equilibrium, but a space-dependent quantity. The output variable is the concentration profile, . Thereby, the LEMC simulation establishes the relation between and necessary to apply the NP equation. The LEMC method correctly computes volume exclusion and electrostatic correlations between ions, so it is beyond the mean-field level of the PNP theory that is routinely applied for nanopores.

The NP equation is solved iteratively with the LEMC simulations in an iteration procedure that ensures that the continuity equation, , is satisfied (more detail is found in the Appendix and in the ESI). The resulting NP+LEMC technique provides a solution for the statistical mechanical problem (e.g., the relation) on the basis of particle simulations, while it still gives an approximate indirect solution for the dynamical problem through the NP equation. Direct simulation of ionic transport in the implicit-water framework is commonly done by the Brownian Dynamics method [20, 21, 22]. The advantage of the NP+LEMC technique is that we can easily handle cases where direct sampling of ions crossing the pore is problematic such as the case of the bipolar nanopore, through which currents are small due to the depletion of ions in the oppositely charged zones.


The questions that we pose in this study are the following. Is the implicit-water model a good approximation to the explicit-water case from the device point of view? Do the two models provide the same device behavior, namely, have the same transfer function for the device? The transfer function describes the relation between the control signal (input) and the response (output) that the model gives to the input signal. Strictly speaking, this is the electrical current as a function of voltage that drives the current. We can define transfer functions, however, in a more general way: the electrical current as a function of experimentally tunable parameters such as pore charge and bulk concentrations.

The transfer functions are the same on the two modeling levels because we study the same device with the two models. The implicit model is justified on the basis of its ability to reproduce the transfer functions given by the detailed model. The MD model, therefore, is the gold standard in this comparison because it contains all the molecular details (explicit water) that are missing from the implicit-water model. If we consider the device as a black box model, we are satisfied with just answering this question. As the dimensions of the devices are getting smaller, however, we want to open the box and look inside to understand how it works on the level of molecules, which is necessary to design new devices.

Even if the explicit- and implicit-water models work similarly on the device-level, opening the box reveals that there are profound differences between the two models due to the presence or absence of explicit water molecules. In this case, new questions arise. How can the transfer functions be similar despite the differences? Specifically, why can water molecules be averaged into an implicit background even in confined spaces? Using a reduced model we may be ignoring crucial details. Looking at too much detail, however, can result in missing the forest for the trees. Unknown errors from inadequately calibrated parameters of detailed models may produce ambiguous artifacts. What are the details that cannot be coarse-grained and should be included in the reduced model to reproduce the device function?

Bipolar nanopore models

To address these questions, we have chosen a bipolar nanopore as a test system because it has a very characteristic transfer function: the pore rectifies electrical current. The output is the current given as a response to the input, the voltage. The current is markedly different for different signs of the voltage due to the asymmetrical surface charge distribution on the pore wall. The sign of the surface charge changes from positive to negative along the central axis of the pore in this study – we call these N and P regions, respectively, according to the nomenclature of the semiconductor literature (Fig. 1). In the N zone (with positive surface charge), the anions are the counter-ions and therefore the majority charge carriers, while the cations are the co-ions and the minority charge carriers. In the P zone, the situation is reversed. Pore regions with opposite surface charges can be achieved by chemical modifications in the case of PET nanopores by transforming carboxyl groups into amino groups by a coupling agent [23], for example. This modification can even be reversible by making the chemical moieties redox-sensitive [24]. The pore charges in the N and P zones do not drive ion tranport; rather, they modulate the concentrations of the various ionic species in the respective zones. The driving force of drift-diffusion of ions is the difference of the concentrations and/or applied voltage on the two sides of the membrane [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 1, 35, 36, 37, 15, 16, 38, 39].

Fabrication of such synthetic nanopores have become possible recently [40, 41, 42, 43, 44, 45, 46, 47]. Their common property (that distinguishes them from micropores) is that their radii are smaller than (or similar in size to) the characteristic screening length of the electrolyte (NaCl in this study) that fills the channel. The dimensionality of these systems makes them central building blocks of devices for various applications from DNA sequencing [48] through switches [49] to chemical and biosensing [50, 51, 52, 53]. Rectification appears (and the nanopore can be called a nanofluidic diode) if the nanopore is asymmetrical in geometry (e.g. conical nanopores) or in charge distribution (e.g. bipolar nanopores studied in this work). What makes the bipolar nanopore more suitable to be the test system of our multiscale modeling study is the fact that the mechanism of rectification is mainly electrostatic in nature and that it rectifies even if the dimensions of the pore are small (1 nm radius and 6 nm length in this work) and therefore suitable to study with MD.

Although bipolar nanopores have been studied with the PNP theory extensively [54, 36, 23, 55, 56, 2, 57, 58, 59, 60, 61, 62, 63, 64], we are not aware of any paper where this system is examined by molecular simulations apart from our previous study for a rectifying bipolar ion channel [65] and our parallel study where we compare to PNP [66]. We are aware only of a few MD studies [67, 68, 69, 70, 71, 72] for other types of nanopores.

In this study we pay special care to modeling the nanopore in the MD and the NP+LEMC calculations as similarly as possible, so that the major difference between the two systems is the treatment of water, compared to which other dissimilarities are minor. This makes this work a case study for studying the differences between the explicit and implicit water models. We connect the two modeling levels by estimating the diffusion coefficients of ions in the pore (for the NP+LEMC calculations) from the MD simulations. We fit the parameters of the functions so that NP+LEMC results for the ionic currents reproduce MD results as closely as possible. In this study, we choose the simplest assumption for the profile: it is a piecewise constant function that is different inside the pore and in the bulk: and . The difference of the and values accounts for all those effects that are different in the pore and in the bulk including different structure of water molecules around the ions and the presence of the pore (confinement and being trapped in electrostatic wells created by the pore charges). We fit the values for one state, fix them, and check whether they are appropriate in other cases too.

We work with a cylindrical nanopore model whose dimensions are relatively small, only nm in radius and nm in length. This radius is quite close to that attainable in experiments at the tip of conical nanopores. The length is far from the experimentally realistic m length scale, but we are restricted by the fact that we use particle simulations. The small pore radius allows the overlapping of the double layers formed at the pore walls at the 0.5 and 1 M concentrations used in this study (wider pores are considered elsewhere [66]). Large concentrations were necessary to achieve a reasonable sampling in the MD simulations. Small pore radius is also useful for the purpose of this study because it amplifies the role of water molecules.

Figure 1: Illustration of the geometry and the rectification mechanism. Top panels: Nanopore structure made by the VMD package [73]. The pore is formed by a carbon nanotube (CNT) between two carbon nanosheets (CNS) defining the membrane. Distance of CNS atoms in the two sheets (width of membrane) is nm, while the distance of CNT atoms from the pore axis (pore radius) is nm. The pore is made charged by placing partial charges on the C atoms of the CNT. The charges are positive on the left hand side (blue dots, N region) and negative on the right hand side (red dots, P region) of the pore. The value of partial charges is chosen so that the average surface charge density corresponds to a prescribed value, . The series of spheres inside the pore represent trajectories of Cl ions through the pore plotted in 10 ps time intervals. The Cl ion was randomly chosen of those that crossed the pore. The trajectories are from simulations for M NaCl and nm at mV voltages (200 mV is ON – left panels, while -200 mV is OFF – right panels). Middle panels: randomly chosen snapshots (from a video available in the ESI) with blue and red spheres representing Na and Cl ions, respectively. The depletion zone of Na is formed in the entire N zone, not only in the PN junction in the middle (the same is true for Cl in the P region). The depletion zones are deeper in the OFF state. Bottom panels: illustrative MD concentration profiles (more detailed profiles are shown in Fig. 3).

The surfaces of the pore and the membrane are smooth. This is advantageous because we can evaluate the results more easily, we can avoid artifacts from surface roughness, and we can create toy models for both the MD simulations (using a carbon nanotube (CNT) as template, see Fig. 1), and the NP+LEMC calculations (using hard walls, see Fig. SI 12). This way, the results of the two methods are comparable and the main difference between them is the treatment of water.

Results and Discussion

Results for the device level

First, we looked at the nanopore from the device perspective, namely, we studied how the response function behaves in the two modeling frameworks. Because MD sampling is weak for small concentrations and voltages, we restrict ourselves to voltages mV and concentration M in the main text. The effect of the bulk concentration and voltage-dependence (current-voltage curves) is found in the ESI showing mainly NP+LEMC results with a few MD data (Figs. SI 1 and SI 5). For a full comparison between the two models, we focus on the dependence of currents on the surface charge density on the nanopore wall, . We place surface charge on the left hand side of the pore wall and on the right hand side (see Fig. 1), therefore, the value of characterizes the strength of the polarity of the pore. We define rectification as the ratio of the magnitudes of the ON (200 mV) and the OFF (-200 mV) currents (ground is on the left hand side).

Figure 2 shows currents as functions of for bulk concentration M (the analogous figure for M is Fig. SI 2 in the ESI). The MD and NP+LEMC currents coincide for /nm in the ON state, because we fitted the ionic diffusion coefficients inside the pore, , so that NP+LEMC results reproduce the MD currents in this case. We obtained the values ms and ms from this adjustment. Then, we fixed these diffusion coefficient values and never adjusted them again; they remain the same for other surface charges, voltages, and concentrations.

Figure 2: Currents as functions of surface charge density for M for the ON (200 mV) and OFF (-200 mV) sign of the voltage. Lines and symbols represent NP+LEMC and MD results, respectively. Blue circles indicate the state point at which the diffusion coefficient values, and , are fitted to the MD currents, therefore, MD and NP+LEMC results coincide in this case.

The relevant questions from the point of view of multiscaling is whether the NP+LEMC calculations using these fixed values can reproduce the MD data for other cases too, namely, in the OFF state, for other surface charges, and other concentrations. Posing the question in a different way; are these values transferable? If they are transferable, the reduced model “knows” the basic physics necessary to describe the device behavior, so no additional fitting is needed. This means that the approximations built into the reduced model smear only the “unimportant” degrees of freedom into response functions ( and ), but leave the “important” degrees of freedom unaffected. A crucial element of multiscaling is identifying the “important” degrees of freedom by doing the calculations at the various resolutions and comparing the results.

Let us examine Fig. 2 in more detail. When the pore is uncharged ( /nm), the currents at the positive and negative voltages are the same (apart from statistical uncertainties; see ESI for details), because the pore is symmetric. When is increased, the currents in the ON and OFF states increasingly deviate. The current (in absolute value) increases in the ON state, while it decreases in the OFF state, i.e., rectification increases with increasing . The MD and NP+LEMC curves are in qualitative agreement regarding this trend. This implies that the profiles as constructed in this work (experimental value in the bulk, , while a constant, fitted, and fixed value in the pore, ) modulate current in a way that agrees with the explicit-water MD model.

Figure 3: Concentration profiles as obtained from MD (left panels) and NP+LEMC (right panels). In the top panels, the Na (solid lines) and Cl (dashed lines) profiles are shown for the ON (black lines) and OFF (red lines) signs of the voltage for M and /nm. In the bottom panels, the ratio of concentration profiles in the ON and OFF states are shown for various values of (see the numbers near the curves). Concentration profiles have been computed in the same way in the two models: the average number of ions in a slab has been divided by the effective volume available for the ions. Inside the pore, the same cross section was used (radius 1 nm) in both cases.

If we want to look into the black box, we need to analyze space-dependent concentration and potential profiles. From the point of view of current, which is the integral of the flux density over the cross section of the pore, the -dependent concentration profiles (averages over cross sections) are the most relevant profiles. They describe which ion is present in which region (N and P) of the nanopore and are the first-order determinants of the current.

Figure 4: (A) Structure of water molecules (turquoise), Na (blue), and Cl (red) ions as obtained from MD simulations in the ON state for /nm and M. Coloring is designed to overemphasize peaks in order to show structures. (B) Radial profiles by averaging over the P region ( nm, top panel) and the N region ( nm, bottom panel) are plotted with symbols connected with lines (same colors as in panel A). Water profiles are divided by 20 to make them comparable to ion profiles. The curves obtained from the implicit-water model are also shown (dashed lines).

The top panels of Fig. 3 show concentration profiles obtained from MD (left panels) and NP+LEMC (right panels) calculations. The two models describe the mechanism of rectification in qualitative agreement. As a primary effect, the pore charges produce depletion zones for the respective ions, e.g., the positive pore charges on the left hand side (N zone) result in depletion zones in the Na profiles (the same is true for the Cl profiles in the P zone). The Na ions, therefore are co-ions in the N zones, while the Cl ions are counter-ions there (reversed in the P zone). In this system, diffusion is limited for both ionic species by their respective depletion zones because these zones are the largest resistance elements if we imagine the ionic pathway as resistors connected in series in an equivalent circuit; if an ion species is not there, it cannot conduct. Depletion zones, therefore, are more important than peaks.

As a secondary effect, the applied field modulates the effect of pore charges. It increases the concentrations (of both ions) in the ON state, while it decreases them in the OFF state. Depletion zones of the ions, therefore, are deeper in the OFF state than in the ON state. This is the basic reason of rectification. It is important to stress, however, that the peaks are also lower in the OFF state. The depletion zones are not independent from the peaks. Co-ions and counter-ions are strongly correlated electrostatically, therefore, co-ions are present in their depletion zones because the counter-ions are present there.

As far as the -dependence is concerned, ON concentrations increase, while OFF concentrations decrease with increasing as shown by the bottom panels of Fig. 3, where the ratio of the ON and OFF concentrations is shown for different values (the concentration profiles themselves are in Fig. SI 4). The ratio increases with increasing . The ratio is larger in the case of the MD data (corresponding to deeper depletion zones), which is reflected in the smaller OFF currents compared to the NP+LEMC results. Better agreement between MD and NP+LEMC could be achieved by adjusting different diffusion constant values in the N and P zones (different mobilities in the N and P zones are implied by the top panels of Fig. 1), but this is beyond the scope of this study.

The results obtained for the operation of the nanopore as a device (and the accompanying profiles) are quite similar on the two modeling levels. This is surprising given that the treatments of water in the two models are so profoundly different. Consequences of the differences can be revealed if we open the black box even wider and dive into details beyond the cross-sectionally averaged profiles.

Results for the molecular level

The presence of water molecules and their effect on ions can be seen by plotting the full profiles. Figure 4A shows contour plots for the ON state as obtained from MD simulations (the plot for the OFF state is seen in Fig. SI 7). A clear layering structure as a function of the coordinate is present. There are two distinct peaks of water, one is near the pore wall solvating the pore charges. The ions’ behavior depends on the region in which we observe them. Cl ions have a large peak just “behind” the solvating water layer in the N region, while they rather accumulate in the pore center in the P region. Na ions, on the other hand, show the opposite behavior, exhibiting two peaks in the P region.

Figure 5: Electrical potential profiles as obtained from MD (left panel) and NP+LEMC (right panel) simulations. The profiles do not include the applied potential. The MD profiles include the contributions of fixed charges (on the pore wall), ions, and water. The MC profiles include only the effect of fixed charges and ions divided by . Results are shown for the ON (top halves) and OFF (bottom halves) signs of the voltage. The profiles have been computed from the average charge profiles provided by the simulations from solving Poisson’s equation (see Appendix and Fig. SI 8).

These effects are better (and more quantitatively) observed in Fig. 4B, where the profiles are plotted averaged over the -coordinate in the N region (top panel) and in the P region (bottom panel). Here, the NP+LEMC profiles for Na and Cl are also shown. Although they are unstructured, they reproduce the relative quantities of the two ion species in the respective regions.

Therefore, despite the differences in the radial structure, the MD and NP+LEMC ionic profiles show basically the same behavior in the dimension. This is the behavior that is relevant for the calculation of the current, and therefore, for the function of the device. The -dependence seems to belong to the class of the “unimportant” degrees of freedom vis-à-vis ion current, because its accurate reproduction is not necessary to properly model the operation of the nanopore as a rectifying diode. The -dependence, on the other hand, is crucially important. This explains why 1-dimensional PNP calculations work so well for nanopores. The water molecules then also seem to belong to the class of “unimportant” degrees of freedom, in the sense that they do not contribute to the current (they have no net charge) and do not affect the axial profiles of ions substantially. Also, the treatment of their dynamic effect on the ions as a friction via the diffusion coefficients in the NP equation seems to be a sufficient approximation.

Their other effect on ions, screening, also works very differently in the implicit and explicit-water models. To perceive these profoundly different ways of screening as given by MD and NP+LEMC simulations, we consider the electrical potential profiles produced by the various species in the two models. The electrostatic potential from water is present in the MD simulations explicitly and provides screening of the potential from ions (and pore charges) in an additive way. In the NP+LEMC simulations, screening is provided by polarization charges that are right on top of the ionic point charges. In effect, this correponds to just damping the electric field of ions by dividing it by (see Eq. 2).

Figure 5 shows potential profiles as contour plots over the coordinates obtained from MD (left panel) and NP+LEMC (right panel) at the same conditions. These potentials do not contain the applied potential, only those produced by free charges in the system: ions, pore charges, and water (if present). Top halves are for the ON state, while bottom halves are for the OFF state. The differences between the two models are striking.

In the implicit-water model (Fig. 5B), only ions can respond to changes in the electric field. The polarization charges move together with the ions so they cannot respond to the electric field independently. The effect of the applied field, therefore, is clearly visible. More ions are available for screening in the ON state, so the potential is smaller (in absolute value) in the ON state. In the OFF state, due to the low concentration of ions, the pore charges dominate the potential.

In the case of the explicit-water model (Fig. 5A), on the other hand, water molecules also produce a counter-field to external effects. The first-order external effects to which both ions and water respond are the applied field and pore charges (these are fixed). Water molecules also respond to the field of ions (and vice versa, of course). Because ions move slowly and they alone (without water) form a low-density (dilute) system, the characteristic speed at which they form the ionic distribution is relatively small. Water molecules, on the other hand, are high-density (liquid) and rotate quickly, so they accomodate themselves to the movement of ions and other external fields fast. The electrolyte, therefore, is a conducting material in which mobile charges are found in high density that respond to external fields and exert a counter-field that produce close to constant potential in the electrolyte.

Figure 6: The two components of the electrical potential of free charges as obtained from MD: potential of the ions and the fixed pore charges (top halves) and potential of water (bottom halves). Results are shown for the ON (panel A) and OFF (panel B) signs of the voltage.

Because the potentials in the implicit-water model are uniformly damped by , the potential varies in a narrow range (). In the explicit-water case, on the other hand, the total potential changes in a much wider range (Fig. 5B vs. A). For this reason, the effect of the relatively small applied field ( mV) is hardly visible.

The effect of the applied field, on the other hand, is better visible on the potentials produced by the ions (together with pore charges) and water molecules separately. These potentials are plotted in Fig. 6. Figures 6A and B show the results for the ON and OFF states, respectively. Water molecules produce an electrical potential (bottom halves) that has an opposite sign compared to the potential of ions and pore charges (top halves).

Differences in the potentials of ions and pore charges in the ON and OFF states are consequences of the differences in the ionic profiles in the ON and OFF states. Water molecules respond easily to the changes in ion distributions (as described above), so their potentials are also different in the ON and OFF states.

Another important difference between the explicit- and implicit-water models is whether ionic double layers are formed near the membrane on the two sides of it. In the implicit-water case they are shown by the splitting cation and anion profiles (Fig. 3 and Fig. SI 9) and by the gradually rising/declining potential profiles (Fig. 5B and Fig. SI 10) on the two sides of the membrane. The signs of the double layer potentials are the opposite in the ON and OFF states and their formation is closely correlated with the relative quantities of counter-ions vs. co-ions in the pore region adjacent to the double layer (see our other papers [65, 66] for details). The presence of the double layers, however, is not necessary to produce rectification. They are absent in the explicit-water MD simulations, still, the different currents in the ON and OFF states are present.


We found definite differences between the results of the explicit- and implicit-water models with respect to (1) how they describe the radial dependence of the particle profiles, (2) whether the double layers on the two sides of the membrane are formed or not, and (3) how the screening of ionic and pore charges by water is done. Still, the -dependence of the ionic concentration profiles is similar in the two models. The effect of the pore charges producing the depletion zones and the effect of the applied field in modulating the effect of pore charges are so powerful that ions respond to them in a similar way no matter whether the water molecules are there explicitly or their screening effect is just a scaling factor in the denominator of the Coulomb potential.

While this is a simple case study, it does answer the paradox of why reduced models work (i.e., reproduce experimental results) in regimes (e.g., confined geometries) where they ought not to work because their molecular-level physics is too approximate. Our results show that some atomic details do not matter for device input-output properties, and that reduced models that reproduce experiments get the overall device physics right, even though they may get some aspects of the molecular-scale physics wrong. That is, reduced models show what device-level physics is necessary/sufficient to predict device characteristics and to design new devices. The results of reduced models should not be overanalyzed, however. The nanoscale-level results they produce are likely to be incorrect in many ways and MD simulations should be used if these details are desired.

For the foreseeable future, however, fully atomistic simulations will probably remain too computationally expensive to regularly predict experiments, and reduced models will remain the only way to compare with experiments in many cases. Our work provides computational insight into how these reduced models should be interpreted.


All-atom MD simulations were performed with the GROMACS (v.5.0.4) program suite [74, 75] using the leap-frog integrator with a 2 fs time step. The CHARMM27 force field [76] was implemented that included 11,000 SPC water molecules, 110-190 Lennard-Jones (LJ) and point charge type anions/cations (depending on the concentration and pore charge density), and LJ atoms for the CNT and CNS. The dimensions of the simulation box were 65.216.8 nm, with periodic boundary conditions (PBC) applied in all spatial directions. The systems were thermostated to 298.15 K by a modified version of the Berendsen (velocity rescaling) algorithm [77]. An ion was considered to cross the channel if it is initially at one side of the membrane and then ends at the opposite side of the membrane after propagating through the channel.

In the reduced model designed for NP+LEMC calculations the membrane and the pore penetrating it are defined by hard walls with which the ions cannot overlap. The cylindrical pore’s radius and length were calculated to mimic the CNT model of the MD simulations as closely as possible on the basis of an estimated distance of closest approach of ions to the carbon atoms. We used the values nm and nm for the pore radius and length, respectively. The fractional point charges have been placed at the same positions as in the CNT model. The Na and Cl ions are charged hard spheres with radii 0.095 and 0.181 nm) immersed in a dielectric continuum () that models the solvent implicitly (Eq. 2). The bulk diffusion coefficients of Na and Cl were ms and ms, respectively. The problem is solved on a discretized grid iteratively; the electrochemical potential is changed until the continuity equation, , is satisfied. The resulting profiles fluctuate around limiting distributions due to statistical uncertainties in the LEMC simulations. The final results are obtained as running averages. Computational time is measured in days (sometimes, weeks) for MD simulations, in hours for the NP+LEMC simulations, but only in seconds for PNP calculations.

Boundary conditions are treated differently in MD and NP+LEMC. In MD, a V/nm homogeneous external electric field was applied along the -direction to achieve mV potential difference. Ions leaving the cell at one end are fed back on the other end due to PBC applied in the dimension. In NP+LEMC, the simulation cell is a finite cylinder due to the rotational symmetry [17, 18, 19]. Two cylindrical compartments on the two sides of the membrane represent the two bulk regions between which the ion transport flows. Both concentrations and electrical potentials are prescribed on the half-cylinders confining these bulk regions.

The potential profiles can be computed either (1) by inserting test charges and sampling the potential on the fly, or (2) by solving Poisson’s equation for the averaged -dependent charge profiles obtained from the simulations (with Dirichlet boundary conditions of assuming zero potential on the closed surface of the confining cylinder). Agremeent between the two methods was excellent for the NP+LEMC data, while it was reasonable for the MD data (see Fig. SI 8). The potential profiles shown in Figs. 5 and 6 have been obtained from the second method.


We gratefully acknowledge the financial support of the Hungarian National Research Fund (OTKA NN113527) in the framework of ERA Chemistry. This material is based upon work supported by the National Science Foundation under Grant No. 1402897 (to D.G.).


  • [1] J. Cervera, B. Schiedt, R. Neumann, S. Mafe, and P. Ramirez. Ionic conduction, rectification, and selectivity in single conical nanopores. J. Chem. Phys., 124(10):104706, 2006.
  • [2] E. B. Kalman, I. Vlassiouk, and Z. S. Siwy. Nanofluidic bipolar transistors. Adv. Mater., 20(2):293–297, 2008.
  • [3] L.-J. Cheng and L. J. Guo. Ionic current rectification, breakdown, and switching in heterogeneous oxide nanofluidic devices. ACS Nano, 3(3):575–584, 2009.
  • [4] Jordan Hoffmann and Dirk Gillespie. Ion correlations in nanofluidic channels: Effects of ion size, valence, and concentration on voltage- and pressure-driven currents. Langmuir, 29(4):1303–1317, 2013.
  • [5] J.-F. Pietschmann, M.-T. Wolfram, M. Burger, C. Trautmann, G. Nguyen, M. Pevarnik, V. Bayer, and Z. Siwy. Rectification properties of conically shaped nanopores: consequences of miniaturization. Phys. Chem. Chem. Phys., 15:16917–16926, 2013.
  • [6] G. Pérez-Mitta, A. Albesa, M. E. Toimil-Molares, C. Trautmann, and O. Azzaroni. The influence of divalent anions on the rectification properties of nanofluidic diodes: Insights from experiments and theoretical simulations. ChemPhysChem, 17(17):2718–2725, 2016.
  • [7] W. Nonner, L. Catacuzzeno, and B. Eisenberg. Binding and selectivity in L-type calcium channels: A mean spherical approximation. Biophys. J., 79(4):1976–1992, 2000.
  • [8] D. Gillespie and D. Boda. The anomalous mole fraction effect in calcium channels: A measure of preferential selectivity. Biophys. J., 95(6):2658–2672, 2008.
  • [9] D. Gillespie. Energetics of Divalent Selectivity in a Calcium Channel: The Ryanodine Receptor Case Study. Biophys. J., 94(4):1169–1184, 2008.
  • [10] D. Boda, M. Valiskó, D. Henderson, B. Eisenberg, D. Gillespie, and W. Nonner. Ion selectivity in L-type calcium channels by electrostatics and hard-core repulsion. J. Gen. Physiol., 133(5):497–509, 2009.
  • [11] D. Gillespie, J. Giri, and M. Fill. Reinterpreting the anomalous mole fraction effect: The Ryanodine receptor case study. Biophys. J., 97(8):2212–2221, 2009.
  • [12] M. Malasics, D. Boda, M. Valiskó, D. Henderson, and D. Gillespie. Simulations of calcium channel block by trivalent ions: Gd competes with permeant ions for the selectivity filter. Biochim. et Biophys. Acta - Biomembranes, 1798(11):2013–2021, 2010.
  • [13] D. Gillespie, L. Xu, and G. Meissner. Selecting ions by size in a calcium channel: The Ryanodine Receptor case study. Biophys. J., 107(10):2263–2273, 2014.
  • [14] D. Boda, É. Csányi, D. Gillespie, and T. Kristóf. Dynamic Monte Carlo simulation of coupled transport through a narrow multiply-occupied pore. J. Phys. Chem. C, 118(1):700–707, 2014.
  • [15] D. Gillespie, D. Boda, Y. He, P. Apel, and Z.S. Siwy. Synthetic nanopores as a test case for ion channel theories: The anomalous mole fraction effect without single filing. Biophys. J., 95(2):609–619, 2008.
  • [16] Y. He, D. Gillespie, D. Boda, I. Vlassiouk, R. S. Eisenberg, and Z. S. Siwy. Tuning transport properties of nanofluidic devices with local charge inversion. JACS, 131(14):5194–5202, 2009.
  • [17] D. Boda and D. Gillespie. Steady state electrodiffusion from the Nernst-Planck equation coupled to Local Equilibrium Monte Carlo simulations. J. Chem. Theor. Comput., 8(3):824–829, 2012.
  • [18] D. Boda, R. Kovács, D. Gillespie, and T. Kristóf. Selective transport through a model calcium channel studied by Local Equilibrium Monte Carlo simulations coupled to the Nernst-Planck equation. J. Mol. Liq., 189:100–112, 2014.
  • [19] D. Boda. Chapter five - monte carlo simulation of electrolyte solutions in biology: In and out of equilibrium. In Ralph A. Wheeler, editor, Ann. Rep. Comp. Chem., volume 10, chapter 5, pages 127–164. Elsevier, 2014.
  • [20] S. H. Chung, T. W. Allen, M. Hoyles, and S. Kuyucak. Permeation Of Ions Across The Potassium Channel: Brownian Dynamics Studies. Biophys. J., 77(5):2517–2533, 1999.
  • [21] W. Im, S. Seefeld, and B. Roux. A Grand Canonical Monte Carlo-Brownian Dynamics Algorithm For Simulating Ion Channels. Biophys. J, 79(2):788–801, 2000.
  • [22] C. Berti, S. Furini, D. Gillespie, D. Boda, R. S. Eisenberg, E. Sangiorgi, and C. Fiegna. A 3-D Brownian Dynamics simulator for the study of ion permeation through membrane pores. J. Chem. Theor. Comput., 10(8):2911–2926, 2014.
  • [23] I. Vlassiouk and Z. S. Siwy. Nanofluidic diode. Nano Lett., 7(3):552–556, 2007.
  • [24] M. Ali, I. Ahmed, P. Ramirez, S. Nasir, S. Mafe, C. M. Niemeyer, and W. Ensinger. A redox-sensitive nanofluidic diode based on nicotinamide-modified asymmetric nanopores. Sensors and Actuators B: Chemical, 240:895–902, 2017.
  • [25] R. M. Lynden-Bell and J. C. Rasaiah. Mobility and solvation of ions in channels. J. Chem. Phys., 105(20):9266–9280, 1996.
  • [26] D. Chen, J. Lear, and B. Eisenberg. Permeation through an open channel: Poisson-Nernst-Planck theory of a synthetic ionic channel. Biophys. J., 72(1):97–116, 1997.
  • [27] W. Nonner and B. Eisenberg. Ion Permeation and Glutamate Residues Linked by Poisson-Nernst-Planck Theory in L-type Calcium Channels. Biophys. J., 75:1287–1305, 1998.
  • [28] M. G. Kurnikova, R. D. Coalson, P. Graf, and A. Nitzan. A lattice relaxation algorithm for three-dimensional Poisson-Nernst-Planck theory with application to ion transport through the Gramicidin Achannel. Biophys. J., 76(2):642–656, 1999.
  • [29] A. E. Cárdenas, R. D. Coalson, and M. G. Kurnikova. Three-dimensional Poisson-Nernst-Planck theory studies: Influence of membrane electrostatics on Gramicidin A channel conductance. Biophys. J., 79(1):80–93, 2000.
  • [30] D. Gillespie, W. Nonner, and R. S. Eisenberg. Coupling Poisson-Nernst-Planck and density functional theory to calculate ion flux. J. Phys.: Cond. Matt., 14(46):12129–12145, 2002.
  • [31] W. Im and B. Roux. Ions and counterions in a biological channel: A molecular dynamics simulation of OmpF porin from Escherichia coli in an explicit membrane with 1 M KCl aqueous salt solution. J. Mol. Biol., 319(5):1177–1197, 2002.
  • [32] D. Gillespie, L. Xu, Y. Wang, and G. Meissner. (De)constructing the ryanodine receptor: Modeling ion permeation and selectivity of the calcium release channel. J. Phys. Chem. B, 109(32):15598–15610, 2005.
  • [33] C. Peter and G. Hummer. Ion transport through membrane-spanning nanopores studied by molecular dynamics simulations and continuum electrostatics calculations. Biophys. J., 89(4):2222–2234, 2005.
  • [34] J. Cervera, B. Schiedt, and P. Ramírez. A Poisson/Nernst-Planck model for ionic transport through synthetic conical nanopores. Europhys. Lett., 71(1):35–41, jul 2005.
  • [35] A. B. Mamonov, M. G. Kurnikova, and R. D. Coalson. Diffusion constant of k inside Gramicidin A: A comparative study of four computational methods. Biophys. Chem., 124(3):268–278, 2006.
  • [36] D. Constantin and Z. S. Siwy. Poisson-Nernst-Planck model of ion current rectification through a nanofluidic diode. Phys. Rev. E, 76(4):041202, 2007.
  • [37] I. Vlassiouk, S. Smirnov, and Z. Siwy. Ionic selectivity of single nanochannels. Nano Lett., 8(7):1978–1985, 2008.
  • [38] C. Song and B. Corry. Testing the applicability of Nernst-Planck theory in ion channels:comparisons with brownian dynamics simulations. PLoS ONE, 6(6):e21204, 2011.
  • [39] J. Cervera, P. Ramírez, S. Mafe, and P. Stroeve. Asymmetric nanopore rectification for ion pumping, electrical power generation, and information processing applications. Electrochim. Acta, 56(12):4504–4511, 2011.
  • [40] Z. Siwy and A. Fulinski. Fabrication of a synthetic nanopore ion pump. Phys. Rev. Lett., 89(19):198103, 2002.
  • [41] Z. Siwy, P. Apel, D. Baur, D. D. Dobrev, Y. E. Korchev, R. Neumann, R. Spohr, C. Trautmann, and K. O. Voss. Preparation of synthetic nanopores with transport properties analogous to biological channels. Surf. Sci., 532:1061–1066, 2003.
  • [42] Z. Siwy, P. Apel, D. Dobrev, R. Neumann, R. Spohr, C. Trautmann, and K. Voss. Ion transport through asymmetric nanopores prepared by ion track etching. Nuclear Instruments & Methods In Phys. Research Section B-beam Interactions With Materials Atoms, 208:143–148, 2003.
  • [43] S. Howorka and Z. Siwy. Nanopores: Generation, engineering, and singly-molecule applications. In P. Hinterdorfer and A. van Oijen, editors, Handbook of Single-Molecule Biophysics, Advances in Chemical Physics, chapter Chapter 11, pages 293–339. Springer, 2009.
  • [44] W. Guo, Y. Tian, and L. Jiang. Asymmetric ion transport through ion-channel-mimetic solid-state nanopores. Acc. Chem. Res., 46(12):2834–2846, 2013.
  • [45] T. Gibb and M. Ayub. Solid-state nanopore fabrication. In Engineered Nanopores for Bioanalytical Applications, pages 121–140. Elsevier BV, 2013.
  • [46] W. Guan, S. X. Li, and M. A. Reed. Voltage gated ion and molecule transport in engineered nanochannels: theory, fabrication and applications. Nanotechnology, 25(12):122001, 2014.
  • [47] Huacheng Zhang, Ye Tian, and Lei Jiang. Fundamental studies and practical applications of bio-inspired smart solid-state nanopores and nanochannels. Nano Today, 8(3):1470–1478, 2016.
  • [48] O. Otto and U. F. Keyser. DNA translocation. In Engineered Nanopores for Bioanalytical Applications, pages 31–58. Elsevier BV, 2013.
  • [49] L. Wang, Y. Feng, Y. Zhou, M. Jia, G. Wang, W. Guo, and L. Jiang. Photo-switchable two-dimensional nanofluidic ionic diodes. Chem. Sci., 8(6):4381–4386, 2017.
  • [50] L. T. Sexton, L. P. Horne, and C. R. Martin. Developing synthetic conical nanopores for biosensing applications. Mol. BioSyst., 3:667–685, 2007.
  • [51] S. Howorka and Z. Siwy. Nanopore analytics: sensing of single molecules. Chem. Soc. Rev., 38(8):2360–2384, 2009.
  • [52] A. Piruska, M. Gong, and J. V. Sweedler. Nanofluidics in chemical analysis. Chem. Soc. Rev., 39:1060–1072, 2010.
  • [53] S. Howorka and Z. S. Siwy. Nanopores as protein sensors. Nat. Biotechnol., 30(6):506–507, jun 2012.
  • [54] H. Daiguji, Y. Oka, , and K. Shirono. Nanofluidic diode and bipolar transistor. Nano Letters, 5(11):2274–2280, 2005.
  • [55] R. Karnik, C. Duan, K. Castelino, H. Daiguji, and A. Majumdar. Rectification of ionic current in a nanofluidic diode. Nano Lett., 7(3):547–551, 2007.
  • [56] I. Vlassiouk, S. Smirnov, and Z. Siwy. Nanofluidic ionic diodes. comparison of analytical and numerical solutions. ACS Nano, 2(8):1589–1602, 2008.
  • [57] R. Yan, W. Liang, R. Fan, and P. Yang. Nanofluidic diodes based on nanotube heterojunctions. Nano Lett., 9(11):3820–3825, 2009.
  • [58] G. Nguyen, I. Vlassiouk, and Z. S Siwy. Comparison of bipolar and unipolar ionic diodes. Nanotech., 21(26):265301, 2010.
  • [59] A. Szymczyk, H. Zhu, and B. Balannec. Ion rejection properties of nanopores with bipolar fixed charge distributions. J. Phys. Chem. B, 114(31):10143–10150, aug 2010.
  • [60] K. P. Singh and M. Kumar. Effect of surface charge density and electro-osmotic flow on ionic current in a bipolar nanopore fluidic diode. J. Appl. Phys., 110(8):084322, 2011.
  • [61] K. P. Singh and M. Kumar. Effect of nanochannel diameter and debye length on ion current rectification in a fluidic bipolar diode. J. Phys. Chem. C, 115(46):22917–22924, 2011.
  • [62] K. P. Singh, K. Kumari, and M. Kumar. Ion current rectification in a fluidic bipolar nanochannel with smooth junction. Appl. Phys. Lett., 99(11):113103, 2011.
  • [63] L. van Oeffelen, W. Van Roy, H. Idrissi, D. Charlier, L. Lagae, and G. Borghs. Ion current rectification, limiting and overlimiting conductances in nanopores. PLOS ONE, 10(5):e0124171, may 2015.
  • [64] M. Tajparast, G. Virdi, and M. I. Glavinović. Spatial profiles of potential, ion concentration and flux in short unipolar and bipolar nanopores. Biochim. Biophys. Acta (BBA) - Biomem., 1848(10, Part A):2138–2153, 2015.
  • [65] Z. Ható, D. Boda, D. Gillepie, J. Vrabec, G. Rutkai, and T. Kristóf. Simulation study of a rectifying bipolar ion channel: detailed model versus reduced model. Cond. Matt. Phys., 19(1):13802, 2016.
  • [66] B. Matejczyk, M. Valiskó, M.-T. Wolfram, J.-F. Pietschmann, and D. Boda. Multiscale modeling of a rectifying bipolar nanopore: Comparing Poisson-Nernst-Planck to Monte Carlo. J. Chem. Phys., 146(12):124125, 2017.
  • [67] A. Aksimentiev, R. K. Brunner, E. Cruz-Chú, J. Comer, and K. Schulten. Modeling transport through synthetic nanopores. IEEE Nanotechnol. Mag., 3(1):20–28, 2009.
  • [68] E. R. Cruz-Chu, T. Ritz, Z. S. Siwy, and K. Schulten. Molecular control of ionic conduction in polymer nanopores. Faraday Discussions, 143:47–62, 2009.
  • [69] E. R. Cruz-Chu, A. Aksimentiev, and K. Schulten. Ionic current rectification through silica nanopores. J. Phys. Chem. C, 113(5):1850–1862, 2009.
  • [70] Q. Chen, L. Meng, Q. Li, D. Wang, W. Guo, Z. Shuai, and L. Jiang. Water transport and purification in nanochannels controlled by asymmetric wettability. Small, 7(15):2225–2231, 2011.
  • [71] T. Gamble, K. Decker, T. S Plett, M. Pevarnik, J.-F. Pietschmann, I. V. Vlassiouk, A. Aksimentiev, and Z. S. Siwy. Rectification of ion current in nanopores depends on the type of monovalent cations – experiments and modeling. J. Phys. Chem. C, 118(18):9809–9819, 2014.
  • [72] Y. Ge, J. Zhu, M. Kang, J. Xian, Q. An, and G. Zhong. Ion current rectification in a confined conical nanopore with high solution concentrations. Mol. Sim., 42(11):942–947, 2016.
  • [73] W. Humphrey, A. Dalke, and K. Schulten. VMD – Visual Molecular Dynamics. J. Mol. Graphics, 14:33–38, 1996.
  • [74] H. J. C. Berendsen, D. van der Spoel, and R. van Drunen. GROMACS: A message-passing parallel molecular dynamics implementation. Comp. Phys. Comm., 91(1–3):43–56, 1995.
  • [75] S. Pronk, Sz. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess, and E. Lindahl. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics, 29(7):845–854, 2013.
  • [76] P. Bjelkmar, P. Larsson, M. A. Cuendet, B. Hess, and E. Lindahl. Implementation of the CHARMM force field in GROMACS: Analysis of protein stability effects from correction maps, virtual interaction sites, and water models. J. Chem. Theor. Comp., 6(2):459–466, 2010.
  • [77] G. Bussi, D. Donadio, and M. Parrinello. Canonical sampling through velocity rescaling. J. Chem. Phys., 126(1):014101, 2007.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description