Ionic Transport in Electrostatic Janus Membranes. An Explicit Solvent Molecular Dynamic Simulation.

Janus, or two-sided, charged membranes offer promise as ionic current rectifiers. In such systems, pores consisting of two regions of opposite charge can be used to generate a current from a gradient in salinity. The efficiency of nanoscale Janus pores increases dramatically as their diameter becomes smaller. However, little is known about the underlying transport processes, particularly under experimentally accessible conditions. In this work, we examine the molecular basis for rectification in Janus nanopores using an applied electric field. Molecular simulations with explicit water and ions are used to examine the structure and dynamics of all molecular species in aqueous electrolyte solutions. For several macroscopic observables, the results of such simulations are consistent with experimental observations on asymmetric membranes. Our analysis reveals a number of previously unknown features, including a pronounced local reorientation of water molecules in the pores, and a segregation of ionic species that had not been anticipated by previously reported continuum analyses of Janus pores. Using these insights, a model is proposed for ionic current rectification in which electric leakage at the pore entrance controls net transport.

Recent efforts to understand the transport of aqueous electrolyte solutions in nanofluidic devices have led to rapid advances. Some of these can be attributed to the availability of improved nanofabrication and characterization techniques, while others have resulted from innovative ideas for energy harvesting and storage.[1,2] The search for devices that maximize ionic selectivity and current rectification while minimizing electric resistance, has led to an array of ingenious designs that are often characterized by symmetry breaking in the axial direction of the pore, and by having radial dimensions that are commensurate with the Debye length of the solution.[3] The so-called ”Janus membranes”, in particular, have garnered considerable attention due to their superior energetic efficiency and the ample design-space they offer.[4−7] Janus membranes consist of two distinct regions with opposite charge. Such membranes exploit the asymmetry created between the two regions to establish unidirectional transport—a feature that is of interest for osmotic energy harvesting in situations where a strong built-in salinity gradient is available (e.g., rivers or oceans). Past reports related to their performance have sometimes been accompanied by theoretical interpretations that rely on solution of the Poisson–Boltzmann (PB) or Poisson–Nernst–Planck (PNP) equations.[5,6,8−10] As helpful as continuum approaches have been for the design of Janus pores, it is unclear whether such models are capable of incorporating all of the relevant physics that arise in nanometer pores.[11,12] Familiar continuum concepts, such as that of a dielectric constant, viscosity, and diffusion coefficient must be revisited in pores, the dimensions of which are comparable to those of the range of interaction between ionic species.

The primary goal of this work is to understand the molecular origin of current rectification in Janus membranes. While this question has been addressed using continuum representations,[5,6] molecular dynamics simulations of rectification have been scarce.[13−15] In particular, all-atom studies have focused on much smaller nanopores than those used in experiments, and have not considered the complex geometries that have been shown to enhance rectification in experiments. The nanopore model considered here is inspired by recent experimental reports, and it exhibits the main features of fully functional nanoscale laboratory devices reported in the literature.[5,6] In the confinement regime considered here (5–10 nm pore diameter), the three-dimensional hydrogen-bond (HB) network of the solution is perturbed but not destroyed. At the same time, electrostatic effects are not too strong because the smallest geometric scale is much larger than the Debye length of the 0.05 M electrolyte solutions in our simulations (1.36 nm). Furthermore, when the concentration is raised to 0.17 M, that Debye length drops to 0.73 nm. This regime is of interest because electrostatic interactions play a key role in the underlying physical processes, but they do not dominate the pairwise interaction energy between molecules, thereby leading to opportunities for more effective control of ionic transport.

Relying on large-scale molecular dynamic simulations, we are able to correctly reproduce many of the experimentally observed membrane characteristics, including the current–voltage response, the charge density distribution, and the ionic current rectification (ICR) factor. Having established the ability of the pore model adopted here to reproduce key experimental features, we are able to identify the underlying molecular processes responsible for the observed behavior.

Our results suggest that the ICR in Janus membranes is linked to the ability of the electric field induced inside the membrane to pull in ions from the reservoir. The asymmetry of the membrane is such that, when a negative bias is applied, the induced field is not able to reach the ions in the reservoir, thereby inhibiting ionic transport. Beyond the asymmetrical ionic distribution inside the pore in response to an applied bias, we find that water also reacts asymmetrically. Water dipoles align strongly with the electric field produced by the charged walls, and also with the externally applied field. Importantly, the dipole alignment is found to be strongly anisotropic, and to depend on the local ionic concentration. The insights provided by the simulations presented here should be of considerable assistance for the improvement of continuum theories used in the description of electrolytes in nanopores, and for the design of devices with better ICR characteristics.

Results and Discussion

We start our discussion with a general characterization of the Janus membrane, and then examine the molecular mechanisms that influence its behavior. Perhaps the most important indicator of performance is the ionic current that is generated in response to an applied potential, shown in Table . Note that a nonlinear response to the applied potential is desirable, and is the reason for the asymmetrical design of the membrane.[1,4,10,16] A qualitative comparison of our results with those reported in the literature[5,6,17] reveals that all the major features observed in experiments on Janus membranes are reproduced by our molecular model, including the ionic selectivity for sodium over chloride, and a significant rectification factor, defined as the ratio between the ionic current of the active state and the inactive state.

Table 1

Ionic Current Generated in the Configurations Considered in This Study, Including Two Concentrations (C) and Three Applied Voltagesa

C [M]	I [nA] activated	I [nA] 0 V	I [nA] inactivated	R (I+/I–)
0.05	3.037	0.022	–0.341	8.91
0.17	6.079	0.048	–1.858	3.27
I enhancement	2.00		5.45

The rectification factor R is calculated by dividing the current obtained for opposite external biases. The current enhancement is defined as the ratio between the current obtained for the two concentrations at the same bias.

The current rectification factor is a crucial indicator of the membrane’s performance. It is also independent of the number of pores, thereby allowing for direct comparison between our results and experimental measurements. Our simulations show a rectification factor of 8.91 at a 0.05 M electrolyte concentration, and one of 3.27 at a 0.17 M concentration, which are consistent with the experimental value of 7.1 reported by Zhang et al. for a similar electrolyte concentration.[5] Other studies on Janus membranes of similar design have reported rectification values of comparable magnitude.[6,17] In agreement with experiments, our simulations show a drop in rectification with increasing electrolyte concentration, as evidenced by the 37% reduction reported in Table when the reservoir concentration is changed from 0.05 M to 0.17 M.

We next consider the average distribution of ions in the absence and presence of an applied electric field. The mechanism for ionic current rectification is generally understood in terms of the ion accumulation and depletion model initially proposed by Woermann for conical pores,[18] and later extended to generic asymmetrical channels.[9,10,16,19,20]Figure left column shows the ionic concentration profiles along the radial and axial coordinates for both the active (applied electric field, 0.0333 V/nm) and inactive versions of the pore (external electric field, −0.0333 V/nm), along with the unbiased equilibrium state. Consistent with the current understanding of current rectification, we find the inactive state to exhibit a reduced ionic content with respect to the active state.[5−7,15] However, in contrast to the Woermann model of ICR, we find that the pore is never fully depleted of ions with respect to the bulk. Instead, the concentration of ions inside the channel is always greater than that in the reservoir.

Figure 1

2D ionic (Na+Cl) concentration distribution (left column) and the radially averaged concentration individual sodium and chloride projected onto x (right column) are shown. The top row is the equilibrium state at 0 V, the middle row is the activated state, while the bottom row is the inactivated case. The 2D ion concentration plot is for a reservoir concentration of 0.05 M. For the higher 0.17 M case please see the Figure S2 in the Supporting Information.

In Figure , the three right-hand column plots show a comparison of the ionic distribution inside the pore for the two electrolyte solution concentrations considered here. For clarity, the ionic concentration has been projected only over the axial coordinate, and split into the anion and cation contributions. The addition of ions does not translate into a substantial change in the distribution of ions within the pore. The bulk of the extra ions is accommodated by the reservoir, raising the total concentration in that region from about 0.05 M to 0.17 M. This implies that the screening of the charged walls of the pores is barely altered. Figure compares the sodium radial concentration profiles in the narrow pore for the two reservoir concentrations. Both configurations exhibit a similar profile, even in the presence of very different Debye electric double layer thicknesses.

Figure 2

Radial concentration distribution of sodium ions in the narrow channel, obtained by averaging over its middle portion 30 nm < x < 50 nm, as a function of r. The results for a reservoir salt concentration of 0.05 and 0.17 M are nearly identical, even though the two Debye lengths, shown using dashed lines, are different.

While most of the pore interior is unaffected by the higher electrolyte concentration, there is a noticeable uptick in the local ion concentrations near the sharp edges at the entrance of the narrow pore, which exposes the wall charge that would have been screened in the equilibrium state. This leads to a significant drop in ϕ around a value of x of 60 nm for the active state when the cations approach from the right, and around 20 nm for the cations moving from left to right in the inactive pore. This can be appreciated in the radially averaged local electric (ϕ) in Figure . It shows a symmetric ϕ profile when adjusted for the direction of ionic movement. These results also confirm that sodium ions are the dominant charge carriers, even if the net charge of the lumen is positive. To show both cation dominance and symmetry more clearly, we also present the radially averaged electric field along the x direction (E) in Figure . The peak location and width are set by the sodium ions for both nonequilibrium states of the channel.

Figure 3

ϕ averaged over the radial coordinate, as a function of x. We show six cases: two reservoir electrolyte concentrations, in the activated, inactivated, and equilibrium states. Also shown is an illustration of the direction taken by the dominant current carrying sodium ions. ϕ has been computed with a reference of 0 at x = 0 and integrating the electric field reported in Figure .

Figure 4

Radially averaged local electric field along the x direction is shown as a function of the axial coordinate x for the active and the inactivate states. At equilibrium, the electric field can be barely detected and is not included here. For reference, the direction of sodium ion transport is also included. The curves shown here are for the reservoir concentration of 0.05 M. For the 0.17 M result, readers are referred to Figure S3. The electric field is calculated from the electrostatic force generated on the ions. This is converted to an electric field via knowledge of the charge. The spatial (r, x) electric field is obtained by averaging over multiple ion instances crossing that location.

There is no asymmetry in the underlying field; the rectification is driven by the local concentration changes that arise from the electrical double layer near the entrance of the narrow section of the pore. Measuring these local concentrations reveals a 2-fold rise in the active case and a 5-fold increase in the inactive state when the reservoir concentration is increased from 0.05 M to 0.17 M. These increases are directly correlated with the higher ionic current reported in Table . To understand ICR further, we also consider the possibility that the asymmetry between the active and inactive modes of the pore might be related to diffusion. To that end, we report the radially averaged axial sodium velocities (V) along x in Figure . We once again recover the symmetry, and the qualitative features of the sodium velocity are consistent with the electric field reported in Figure , underscoring that enhanced diffusion does not account for the magnitude of R of nearly 9 reported in Table .

Figure 5

Sodium axial velocity, averaged over the radial coordinate, as a function of position x for both the active and inactive states of the channel. The ions in the reservoir and in contact with the outer wall are excluded from the average. A complete 2D velocity snapshot is presented in Figure S4.

Our proposed electric-field leakage model for ICR is based on both the electric force and the concomitant acceleration that starts well before the ion-depleted regions of the small section of the pore. In the active pore, ions are brought to the entrance from as far as 5 nm away, in the reservoir, by the electric field generated by the pore when voltage is applied. Similarly, in the inactive pore, ions are accelerated from remote positions but, this time, in the interior of the positively charged chamber, where the sodium concentration is much lower (even lower than in the bulk). Thus, the impact of ionic concentration on current is through an increase in the probability of crossing a barrier, and not necessarily a change in the barrier height.

To better illustrate the dynamics of the Janus membrane operation we show in Figure the path taken by a sodium ion as it traverses the nanopore. We can observe a Brownian motion mechanism, aided by the applied electric field, through which the sodium ion approaches the membrane entrance. There, it spends a considerable amount of time before it crosses over to the negatively charged section of the pore. Once in that section, it is rapidly shuttled to the wall and dragged toward the wider section, where it is pushed to the center of the channel.

Figure 6

Representative trajectory of a sodium ion, the dominant charge carrier, across the membrane. It begins with classical Brownian motion in the reservoir. As it enters the narrow pore, the ion is driven to the wall and dragged along its surface. As the ion approaches the wider portion of the pore, it is forced to the center of the channel. The return to Brownian motion on the other side of the pore is not shown for conciseness.

We now turn our focus toward the explicit role of water in a system with dimensions comparable to those used in laboratory experiments. It is well-known that water behavior deviates from that of the bulk under extreme confinement conditions. In the context of nanofluidics, this level of confinement is often understood to occur below one nanometer, where the size of the system is comparable to molecular dimensions.[11,21,22] However, the size of the pore considered here is 1 order of magnitude above that scale, and therefore in a regime that has been assumed in the past to approach the continuum limit. In that regime, classical theories are expected to be valid, but have been less explored by MD simulations.[12] In what follows, we only highlight some of the most salient features of the structure and dynamics of water in our simulated Janus membrane.

One of the main reasons that water behavior is difficult to predict is that molecular interactions propagate well beyond the size of a single molecule through a complex hydrogen-bonding network. Furthermore, many of the unusual phenomena observed in nanofluidic devices occur in pore sizes that are too small to accommodate a fully three-dimensional HB network.[23−25] Our results indicate that the number of hydrogen bonds per molecule inside the pore is similar to the value computed in the reservoir.

The average dipole moment under negative and positive applied potentials is shown in Figure . One can appreciate a clear tendency for water molecules to align, with both the external electric field and the electric field emanating from the charged walls. The anisotropic nature of the order spans a length that is many times larger than the size of a water molecule. There is some precedent for the dipole behavior found here in the literature.[26−29] Past studies also proposed a reduction of the dielectric constant in the direction perpendicular to the confinement, and an increase in the parallel direction. The strength of the electrostatic interactions and, therefore, the dynamics of ionic motion, could be significantly affected by the distinctive orientation of water at the entrance of the smaller section of the pore.[30] Importantly, although the dipole moment is altered in these regions, the average number of HB per molecule remains at 3.47, which almost exactly matches the average all across the system, including the reservoir, of 3.44.

Figure 7

A 2D slice showing the dipole moment of water when the reservoir concentration is 0.05 M for both the active (a) and inactive (b) pore. There is strong alignment next to the charged wall in both cases, but the peak along x occurs at the respective depletion regions. Results for the 0.17 M case are shown in Figure S5.

The results in Figure allow us to make two additional, important observations. First, as seen in panel b, the orientation of water molecules starts inside the pore, but extends well into the reservoir, for about 5 nm, highlighting the need for a reservoir when this kind of system is modeled either by molecular simulations or continuum approaches.[31] Second, the zones displaying the strongest water alignment in Figure , naturally correlate with the regions of augmented electric field displayed in Figure , which are also the depletion regions seen in Figure .

For completeness, in Figure , we show how the water velocity developed in the radial direction is related to the sodium velocity that drives it. From this figure, one can appreciate a clear correlation between the two velocities. In fact, away from the depletion region, water velocity exactly matches the sodium velocity. Figure b shows how the proportionality is maintained even at different electrolyte concentrations and applied biases, for the region (30 nm < x < 50 nm). These data are consistent with previous reports for pores of smaller sizes, which have noted the concurrence of ionic current and the underlying electro-osmotic flow.[15,32,33] This behavior can be explained by the considerable influence that one of the ionic species (sodium in our case) has over the transport across the channel.

Figure 8

Water velocity averaged along the radial direction and projection onto the axis of the pore, along with the sodium result for reference, as a function of x in panel a. The full lines are sodium velocities: blue lines are active 0.05 and 0.017 M (darker), the red lines are inactive 0.05 M and 0.17 M (darker). The dashed lines are water velocities for inactive (pink, both concentrations) and light blue for active (both concentrations) pores. For ease of comparison with Figure , the velocities computed between 30 nm < x < 50 nm for both sodium and water are shown in panel b. There is no lag in the sodium motion relative to the underlying solvent within the narrow pore. The complete 2D water velocity information for a reservoir concentration of 0.05 M is given in Figure S6.

Beyond the dynamics (via coupling of sodium and solvent velocity) and complex electrostatic effects (e.g., the nontrivial dipole alignment), there is another way in which the often overlooked presence of water is manifest in all the interactions taking place within the system. Figure shows the immediate vicinity of the charged wall. The sodium ions belonging to the so-called “stern layer” or “stagnant layer” are not in direct contact with the charged particles of the walls of the pore, but instead are located precisely between the first and the second layers of the water hydrating the surface. In addition to the “static” ion distribution, the nonzero water velocity that arises from the surface roughness drags the ions in the stern layer along with it, thus impacting the dynamics of the system. The fact that water mediates interfacial interactions certainly introduces more complexity to the already intricate structure and dynamics of the electrical double layer.

Figure 9

Water mediating the contact between the ions and the walls of the pore. The radial concentration profile for water and ions in the inset reveals that sodium ions sit precisely between the first and the second layer of the water’s oxygen atoms. These results are for the 0.05 M reservoir and active state simulation.

Conclusions

The results of equilibrium and nonequilibrium molecular dynamics simulations of a Janus pore in the active and inactive states are in agreement with a number of experimentally determined observables from asymmetric membranes, serving to underscore the validity of the models adopted here.

The electric current and rectification reported in this work can be explained in terms of an electric-field leakage model. The “leaking” electric field arises from the applied electric field, and rearranges ions to expose wall charges. The crossing of these depletion regions at the entrance of the narrow pore, along the direction of the dominant (current carrying) sodium ions, is proposed to occur through the induced electric field. This electric field leaks into the reservoir for the active case, and pulls in the cations to increase conduction. In the inactive case, the leakage of the electric field does not reach into the reservoir, and instead ends in the positively charged chamber, where the sodium concentration is very small, resulting in a low output. The addition of electrolytes to the bulk (for both directions of applied electric field) should increase sodium availability, thereby enhancing electric throughput. More specifically, the enhancement depends on the increase of the local concentration in the narrow pore entrance, which depends on the spatial evolution of the double layer. Hence, to increase ICR, measures could be taken to limit sodium concentration at the narrow/wide pore junction (at x = 10 nm) such as extending the wider section of the membrane, or increasing its positive wall charge density.

Our comprehensive characterization of the water structure and its influence over the membranes’ performance indicates that water dipoles align not only with the externally applied electric field, but also with the field that originates from the charged walls. More importantly, water orientation is found to be highly anisotropic, with a larger average dipole moment that coincides with the ionic depletion regions. Water molecules not only align their dipoles toward the charged walls, but also in a way in which dipole-induced interactions dominate the physics at the surface of the pore. Ions, in contrast, occupy the space between the first and the second layers of water, instead of being in direct contact with the charged particles of the walls.

The results of this work confirm that ICR is not limited to pore radii below the Debye length. We find that the radial ion distribution in the central portion of the membrane is not too sensitive to the electrolyte concentration in the reservoir.

Model and Methods

In the experiments by Zhang et al.,[5] the active part of the Janus pore consisted of two contiguous cylindrical sections. The narrow section had a surface charge of −0.08 C/m2, while the wider section had a surface charge of +0.24 C/m2. We retained what was designed in our simulations but rescaled the length of the pore to 10 nm at the wider section and 40 nm at the narrower part as to preserve the original length and diameter ratio between the oppositely charged sections and still be numerically feasible. In our model a “reservoir” with dimensions 20 nm × 28.5 nm × 28.5 nm is connected to both sides of the pore. A schematic representation of the system is shown in Figure . This Janus membrane is operated with a NaCl 1:1 electrolyte solution. The molecular simulations were conducted under the assumption that the transport properties in the real membrane consist of the sum of the transport of many noninteracting independent pores working in parallel. Under this supposition, the behavior of the membrane is expected to be recovered from the simulation of a single pore. This is justified by the negligible change in conductance for different pore arrangement densities at ionic concentrations similar to ours in the experiments by Zhu et al.[17]

Figure 10

A 2D slice of the model Janus membrane used in our molecular dynamics simulations. The gray spheres are the fixed walls; pink represents water molecules, blue is sodium ions, and green is chloride ions. Only the horizontal walls are charged.

For simulations, the nanopore was carved from an FCC structure with a lattice constant of 0.407 nm. The particles of the walls are modeled with Lennard-Jones parameters σ = 0.32 nm and ε = 0.65 kJ/mol. The small section of the channel is composed of 41 691 particles, each carrying a charge of Q = −0.01504 au to achieve the desired surface charge. Similarly, the large section of the channel contains 17 921 particles, each carrying a charge of 0.04464 in atomic units. The channel is embedded in a box with dimensions provided in Figure , and is filled with approximately 705 000 water molecules. Using pre-equilibrated NaCl solutions, two simulations are performed, with 2234 and 4393 ions, in an overall electroneutral simulation box. For the electrolyte solution, we use the Dang parameters[34] for the ions and SPC/E[35] for the water.

All the MD simulations were performed using GROMACS[36] in the NVT ensemble with periodic boundary conditions. Long-range electrostatic forces were calculated using the particle-mesh Ewald method, with a grid density of 0.16 nm and a real space cutoff at 1.2 nm. The temperature of all simulations was kept constant at 300 K using a velocity-rescaling algorithm coupled every 0.1 ps with no center of mass motion removal. The positions of all the particles of the walls were restrained using harmonic potentials of 9000 kJ·mol–1·nm–2 over the entire course of the simulations. A uniform electrostatic field of (±)0.0333 V/nm, corresponding to active and inactive, states, respectively, was applied in the axial direction to produce the ionic current.

The entire system was first equilibrated without external bias for a period of 30 ns. After that, a production simulation was run for another 30 ns to generate the reference equilibrium state of the pore. Starting from the equilibrium state, intermediate simulations were conducted at each different potential for a period of 30 ns until steady state was reached. The simulations were continued for an additional 30 ns period for the final data collection. The steady state of the final runs was determined by inspecting the ionic current evolution over the course of the simulations (see Figure S1 ).